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Abstract 

This  dissertation  demonstrates  how  the  symmetric  convolution-multiplication  property  of 
discrete  trigonometric  transforms  can  be  applied  to  traditional  problems  in  image  reconstruction 
with  slightly  better  performance  than  Fourier  techniques  and  increased  savings  in  computational 
complexity  for  symmetric  point  spread  functions.  The  fact  that  the  discrete  Fourier  transform 
diagonalizes  a  circulant  matrix  provides  an  alternate  way  to  derive  the  symmetric  convolution- 
multiplication  property  for  discrete  trigonometric  transforms.  Derived  in  this  manner,  the  sym¬ 
metric  convolution-multiplication  property  extends  easily  to  multiple  dimensions  and  generalizes 
to  multidimensional  asymmetric  sequences.  The  symmetric  convolution-multiplication  property 
allows  for  linear  filtering  of  degraded  images  via  point-by-point  multiplication  in  the  transform 
domain  of  trigonometric  transforms.  Specifically  in  the  transform  domain  of  a  type-II  discrete 
cosine  transform,  there  is  an  asymptotically  optimum  energy  compaction  about  the  low- 
frequency  indices  of  highly  correlated  images  which  has  advantages  in  reconstructing  images 
with  high-frequency  noise.  The  symmetric  convolution-multiplication  property  allows  for  well- 
approximated  scalar  representations  in  the  trigonometric  transform  domain  for  linear  reconstruc¬ 
tion  filters  such  as  the  Wiener  filter.  An  analysis  of  the  scalar  Wiener  filter's  improved  mean- 
squared  error  performance  in  the  trigonometric  transform  domain  is  given. 
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TRIGONOMETRIC  TRANSFORMS  FOR 


IMAGE  RECONSTRUCTION 

I.  Introduction 

The  results  of  the  research  presented  in  this  dissertation  demonstrate  new  forms  of  linear 
image  reconstruction  filters  which  employ  a  recently-developed  property  of  trigonometric  trans¬ 
forms.  Image  reconstruction  is  the  process  of  restoring  degraded  images  [40].  An  imaging  sys¬ 
tem  typically  measures  a  blurred,  noisy  image  of  an  object.  The  system  must  then  apply  image 
reconstruction  techniques  to  recover  an  estimate  of  the  object  from  its  blurred,  noisy  version. 
Linear  image  reconstruction  techniques  use  linear  shift-invariant  filters  to  recover  the  ob¬ 
ject  [27].  A  linear  shift-invariant  filter  is  a  system  whose  response  to  an  arbitrary  input  is  com¬ 
pletely  characterized  by  its  response  to  a  single  impulse  [37]. 

The  Air  Force  is  interested  in  image  reconstruction  because  it  has  a  need  to  image  space- 
borne  objects  from  the  ground  [13].  This  problem  is  complicated  by  the  turbulent  atmosphere 
which  has  a  severe  degrading  effect  on  the  quality  of  images  [40]  and  by  the  noise  present  in  im¬ 
age  detection  systems  [17].  Many  existing  nonlinear  image  reconstruction  techniques  are  itera¬ 
tive  [30]  -  [32],  and  require  large  amounts  of  computer  processing  time  with  a  high  degree  of 
computational  complexity.  Existing  linear  methods  [15],  [21],  [27],  [38],  [44],  require  fewer 
computations  than  iterative  methods,  but  suffer  from  poorer  performance.  An  image  reconstruc¬ 
tion  technique  which  has  better  performance  than  existing  linear  techniques,  yet  still  offers  the 
computational  advantage  of  linearity  would  thus  be  of  great  interest  to  the  Air  Force. 
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The  general  purpose  of  this  dissertation  is  to  investigate  the  application  of  trigonometric 
transforms  to  the  two-dimensional  image  reconstruction  problem.  Two-dimensional  trigonomet¬ 
ric  transforms  [26],  [27]  convert  the  pixels  of  a  discretely-sampled  image  into  a  two-dimensional 
series  of  coefficients  which  represents  the  amount  of  information  contained  in  the  image  at  sinu¬ 
soids  of  different  spatial  frequencies.  A  spatial  frequency  is  the  two-dimensional  equivalent  of  a 
one-dimensional  temporal  frequency.  Specifically,  this  research  shows  how  the  symmetric  con- 
volution-multiplication  property  of  discrete  trigonometric  transforms  can  improve  linear  image 
reconstruction  filters. 

Martucci  [34]  recently  developed  the  symmetric  convolution-multiplication  property  for 
the  family  of  discrete  trigonometric  transforms.  The  discrete  trigonometric  transform  family 
consists  of  sixteen  different  one-dimensional  transforms  which  are  even  and  odd-length  versions 
of  type  I  -  rv  discrete  sine  and  cosine  transforms  [26],  [27],  [39].  A  sine  or  cosine  transform  per¬ 
forms  a  similar  operation  on  an  image  as  a  Fourier  transform.  Consider  the  sequence 

=  7{^(«„K2)}  which  is  the  two-dimensional  trigonometric  transform  of  the  sequence 
0(«i,«2)-  The  subscript  T' denotes  that  the  sequence  i9j.(A:,,^2)  exists  in  the  trigonometric 
transform  domain.  The  operator  7{»}  is  a  two-dimensional  discrete  trigonometric  transform 
based  on  any  of  the  sixteen  one-dimensional  discrete  trigonometric  transforms.  The  indices  «[ 
and  «2  correspond  to  an  ordering  of  the  pixels  of  the  image.  The  indices  and  k2  represent  the 
locations  of  sample  points  in  trigonometric  transform  domain  space  which  correspond  to  differ¬ 
ent  spatial  frequencies.  The  values  of  the  sequence  &j(k^,k2)  at  each  point  represent  the 
amount  of  information  contained  in  the  image  which  agrees  with  sinusoids  having  spatial  fi'e- 
quencies  corresponding  to  A:,  and  ^2- 
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Similarly  the  sequence  &p{k^,k{)  =  ?{^(«i,n2)}  the  Fourier  domain  represents  the 

amount  of  information  contained  in  an  image  which  agrees  with  complex  exponentials  of  differ¬ 
ent  spatial  frequencies.  Here  the  subscript  'F'  indicates  the  Fourier  domain,  and  the  operator 
?{•}  represents  the  discrete  Fourier  transform.  The  sequences  (^1,^:2)  and  &p{k^,k2)  are 
both  transform-domain  spatial  frequency  representations  of  the  sequence  0{n^,n2).  One  notable 
difference  between  the  sequences  ^nd  9p{k^,k2)  is  that  if  6(n^,n2)  is  real-valued, 

then  the  sequence  9p{k^,k2)  will  also  be  real-valued,  but  the  sequence  i9^(^,,^2)  will  be  com¬ 
plex-valued. 

There  are  sixteen  different  one-dimensional  discrete  trigonometric  transforms.  Different 
transforms  exist  for  sines  and  cosines,  for  even  and  odd-length  sequences,  and  for  different  types 
distinguished  as  types  I  -  IV  [26],  [27],  [39].  The  four  types  I  -  IV  impose  half-sample  shifts  in 
the  input  or  output  indices  of  the  sine  and  cosine  transforms.  A  type-I  trigonometric  transform 
imposes  no  shift  to  either  the  input  or  the  output  sequence  indices.  A  type-II  transform  imposes  a 
half-sample  shift  to  just  the  input  index.  A  type-III  transform  imposes  a  half-sample  shift  to  just 
the  output  index.  A  type-fV  transform  imposes  a  half-sample  shift  to  both  the  input  and  the  out¬ 
put  indices.  The  fact  that  these  four  transforms  require  either  no  shift  or  a  half-sample  shift  is 
related  to  the  idea  of  the  point  of  symmetry  in  the  symmetric  extension  of  a  finite  sequence.  The 
point  of  symmetry  can  be  either  the  end  point  in  the  sequence  or  a  point  which  lies  one-half 
sample  beyond  the  end  point  of  the  sequence. 

The  discrete  trigonometric  transform  family  lies  at  the  heart  of  many  image  transform 
coding  applications  [27],  [33],  [39].  The  objective  of  image  transform  coding,  which  is  different 
than  the  goal  of  this  research,  is  to  reduce  the  information  content  of  an  image  for  storage  or 
transmission.  Although  very  useful  for  transform  coding,  the  discrete  trigonometric  transforms 
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have  not  proved  very  useful  in  a  wide  variety  of  image  filtering  applications,  because  until  re¬ 
cently  no  convolution-multiplication  property  existed  for  the  entire  family.  The  advantage  of  a 
transform  that  possesses  a  convolution-multiplication  property  is  that  the  property  allows  a  filter 
to  be  implemented  more  efficiently  in  the  transform  domain  through  point-wise  multiplication. 
To  implement  the  effects  of  a  filter  through  convolution  directly  in  the  spatial  domain  requires  a 
large  number  of  shifts,  adds,  and  multiplies. 

The  circular  convolution-multiplication  property  of  discrete  Fourier  transforms  is  well- 
known  [37].  A  symmetric  convolution-multiplication  property  for  the  entire  family  of  discrete 
trigonometric  transforms  has  only  existed  since  1994  [34].  An  earlier  attempt  to  define  a  convo¬ 
lution-multiplication  property  for  the  discrete  cosine  transform  [6]  produced  a  result  limited  to 
the  convolution  of  only  certain  types  of  symmetric  sequences.  It  was  based  solely  on  the  type-II 
even-length  discrete  cosine  transform. 

Martucci  [34]  defines  symmetric  convolution  as  the  form  of  convolution  for  the  entire 
family  of  discrete  trigonometric  transforms  similar  to  the  manner  in  which  circular  convolu¬ 
tion  [37]  is  the  form  of  convolution  for  the  discrete  Fourier  transform.  Circular  convolution  in¬ 
volves  the  point-by-point  shifting  and  adding  of  periodic  replicas  of  the  finite  sequences  being 
convolved.  Symmetric  convolution  involves  the  point-by-point  shifting  and  adding  of  symmetric 
extensions  of  the  sequences  being  eonvolved.  The  symmetric  convolution-multiplication  prop¬ 
erty  states  that  an  inverse  trigonometric  transform  of  the  product  of  the  trigonometric  transforms 
of  two  sequences  yields  the  same  result  as  the  symmetric  convolution  of  the  two  sequences  [34]. 
Mathematically,  the  symmetric  convolution-multiplication  property  states  that  the  sequence 

which  is  calculated  as  c?(«,,«2)  =  '73"’|7i{^(ni,«2)}''72{/?(«i,«2)}}  is  exactly  the  same 
sequence  which  results  fi’om  computing  the  symmetric  convolution  of  the  sequences  ^(«,,«2) 
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and  ;z(«„«2)  directly.  The  operators  72{.},  and  73{»}  are  different  two-dimensional  dis¬ 
crete  trigonometric  transforms.  The  circular  convolution-multiplication  property  for  discrete 
Fourier  transforms,  =  7‘'{7{^(n„«2)}  •  7{Kn^,Th)}],  is  veiy  similar  to  the  symmetric 

convolution-multiplication  property  for  discrete  trigonometric  transforms.  The  operator  7{.} 

again  represents  the  discrete  Fourier  transform.  The  symmetric  convolution-multiplication  prop¬ 
erty  exists  for  forty  different  one-dimensional  cases  of  various  combinations  of  the  sixteen  trans¬ 
forms  in  the  discrete  trigonometric  transform  family.  Symmetric  convolution  is  limited  to  se¬ 
quences  having  underlying  symmetiy.  The  forty  different  cases  cover  many  possibilities  for  both 
symmetric  and  antisymmetric  sequences  which  allows  an  asymmetric  sequence  to  be  decom¬ 
posed  into  its  symmetric  and  antisymmetric  parts  prior  to  being  convolved. 

Previous  applications  of  the  discrete  cosine  transform  to  linear  image  reconstruction  [26], 
[27]  provided  very  good  diagonal  approximations  for  certain  types  of  correlation  matrices.  Many 
linear  image  processing  techniques  require  knowledge  of  how  the  individual  pixels  of  an  image 
are  correlated  with  each  other.  If  the  elements  of  the  vector  9  represent  the  pixels  of  the  image 

^(«,,«2),  the  correlation  between  pixels  is  expressed  by  the  matrix  Rgg  =  66^ .  The  overbar,  ~ , 
indicates  the  expected  value  operator,  and  the  superscript  T' denotes  the  transpose  of  the  vector 
9.  If  the  size  of  the  image  is  iV,  x  7\^2>  then  9  will  be  an  x  1  vector  and  Rgg  will  be  an 
NyN2  X  N^N2  matrix.  The  type-II  discrete  cosine  transform  has  the  property  that  it  very  nearly 
diagonalizes  the  correlation  matrix  of  an  image  whose  pixels  are  highly  correlated  [27].  An  im¬ 
age  with  highly-correlated  pixels  has  few  abrupt  transitions  within  a  small  neighborhood  of  pix¬ 
els  in  the  image.  The  approximate  diagonal  form  of  the  correlation  matrix  for  a  highly-correlated 

image  is  Rg^g^  —  C^e,N^N^^ee^ue,N^Nl•  The  N^N2  x  A^iA^2  matrix  is  a  type-II  discrete 
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cosine  transform  matrix  for  even-length  two-dimensional  sequences.  The  discrete  cosine  trans¬ 
form  of  the  vector  6  is  the  vector  &j.=  The  off-diagonal  elements  of  the  matrix 

approximately  zero.  As  the  pixels  in  an  image  become  more  highly  correlated,  the 
off-diagonal  elements  in  its  transform-domain  correlation  matrix,  approach  zero  [39]. 

One  of  the  results  of  this  research  demonstrates  another  diagonalizing  property  of  discrete 
trigonometric  transforms.  Symmetric  convolution  has  an  equivalent  matrix  representation  as 
d  =  Hsq9.  The  matrix  incorporates  all  of  the  symmetric  extensions,  shifts,  additions,  and 
multiplications  of  symmetric  convolution.  The  resulting  vector  </ contains  the  same  elements  as 
the  elements  in  <i(«j,«2)  which  resulted  from  the  symmetric  convolution  of  the  two-dimensional 
sequences  6{n^,n2)  and  h(n^,n2).  In  the  vector-matrix  representation,  the  sequences  9{n^,n2) 
and  /i(«,,n2)  are  represented  by  the  vectors  9  and  A. 

Trigonometric  transforms  diagonalize  symmetric  convolution  matrices  because  symmetric 
convolution  in  the  spatial  domain  is  equivalent  to  point-wise  multiplication  in  the  transform  do¬ 
main.  The  results  of  this  research  reveal  diagonalizing  forms  for  all  forty  cases  of  symmetric 
convolution.  In  general,  d  =  where  the  matrix  is  diagonal  with  the  vector  T^h 

along  its  diagonal.  This  expression  is  exactly  equivalent  to  d  =  Tf^[Ti9  ©r2^].  where  ©  rep¬ 
resents  a  point-wise  or  Hadamard  product  [22].  The  three  matrices  7j,  T2,  and  are  matrix 
representations  of  two-dimensional  trigonometric  transforms.  The  fact  that  the  above  matrix 
equation  produces  exactly  the  same  result  as  «:/(«, ,«2 )  =  7i'  |7]  {^(«i  ,n2 )} '  {/*(«! ,«2 )}}  is  one 

of  the  main  contributions  of  this  research  [9]. 

The  matrix  form  of  the  s)mmetric  convolution-multiplication  property  allows  for  a  more 
natural  extension  of  the  property  to  asymmetric  multidimensional  sequences  than  the  operator 
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form  of  the  property.  The  ability  to  apply  the  property  to  asymmetric  two-dimensional  sequences 
is  important  because  they  are  the  most  general  class  of  sequences  encountered  in  image  recon¬ 
struction.  The  alternate  derivation  of  the  matrix  form  of  the  property  depends  on  the  fact  that  a 
discrete  Fourier  transform  matrix  diagonalizes  a  circulant  matrix  representing  circular  convolu¬ 
tion  [23]. 

The  diagonalizing  results  for  symmetric  convolution  matrices  presented  here  are  related  to 
the  work  of  Sanchez  et  al.  Her  team  shows  that  diagonalizing  forms  exist  for  the  eight  discrete 
cosine  transforms  [41]  and  the  eight  discrete  sine  transforms  [42].  Their  method  is  to  determine 
the  class  of  matrices  whose  eigenvectors  are  sine  and  cosine  transforms,  but  their  results  reveal 
diagonalizing  forms  for  only  sixteen  of  the  forty  cases  of  symmetric  convolution  ~  one  case  of 
symmetric  convolution  for  each  of  the  sixteen  transforms.  Their  results  therefore  do  not  extend 
to  the  general  case  of  convolving  an  asymmetric  sequence  which  requires  different  convolution 
relations  for  each  part  of  the  sequence's  underlying  symmetry.  The  approach  taken  here  is  to  re¬ 
late  the  discrete  trigonometric  transform  to  the  discrete  Fourier  transform  which  already  pos¬ 
sesses  a  diagonalizing  form  [23].  These  newly  derived  diagonalizing  matrix  forms  then  allow  for 
a  natural  extension  of  all  the  cases  of  the  symmetric  convolution-multiplication  property  to  mul¬ 
tidimensional  asymmetric  sequences. 

The  notion  that  a  type-II  discrete  cosine  transform  matrix  approximately  diagonalizes  the 
correlation  matrix  of  a  highly-correlated  image  and  the  property  of  the  trigonometric  transform  to 
exactly  diagonalize  a  symmetric  convolution  matrix  both  play  important  roles  in  applying  trigo¬ 
nometric  transforms  to  image  reconstruction  problems.  Many  existing  linear  image  reconstruc¬ 
tion  techniques  rely  on  knowledge  of  both  the  correlation  of  the  pixels  in  an  object  being  imaged 
and  the  impulse  response  of  the  system  which  degrades  the  object.  Optical  scientists  and  engi¬ 
neers  refer  to  the  two-dimensional  impulse  response  of  the  degrading  system  as  its  point  spread 
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function.  The  name  arises  from  the  blurring  or  spreading  of  individual  points  comprising  the 
object  [16]. 

As  mentioned  previously,  it  is  computationally  more  efficient  to  calculate  convolutional 
results  in  the  transform  domain.  The  calculations  become  even  more  efficient  if  all  the  transform 
domain  matrices  involved  are  either  diagonal  or  well-approximated  by  their  diagonal  elements. 
Transform  domain  implementations  of  filters  which  rely  on  just  the  diagonal  elements  of  matri¬ 
ces  are  called  scalar  filters. 

Traditional  Fourier  image  reconstruction  filters  have  scalar  forms.  For  an  image  recon¬ 
struction  filter  to  have  a  scalar  form,  both  the  matrix  representing  convolution  of  the  degrading 
point  spread  function  and  the  correlation  matrix  of  the  object  must  have  diagonal  or  approxi¬ 
mately  diagonal  forms.  Discrete  Fourier  transform  matrices  diagonalize  circulant  convolution 
matrices  [23].  Discrete  Fourier  transform  matrices  also  provide  approximate  diagonal  forms  of 
object  correlation  matrices  in  the  transform  domain  [48].  For  a  highly  correlated  image,  the  dis¬ 
crete  Fourier  transform  will  not,  however,  condense  as  much  of  the  content  in  all  the  transform 
domain  terms  onto  the  diagonal  as  a  type-II  discrete  cosine  transform  [20],  [29].  The  lower  en¬ 
ergy  in  the  off-diagonal  elements  in  the  cosine  transform  domain  leads  to  a  better  diagonal  ap¬ 
proximation  than  the  discrete  Fourier  transform. 

The  results  of  this  research  show  for  the  first  time  that  trigonometric  transform  represen¬ 
tations  of  image  reconstruction  filters  also  have  scalar  forms.  The  fact  that  a  type-II  discrete  co¬ 
sine  transform  matrix  approximately  diagonalizes  the  correlation  matrix  of  a  highly-correlated 
image  is  well-known  [27].  More  recently,  Martucci's  symmetric  convolution-multiplication 
property  for  discrete  trigonometric  transforms  [34]  allows  the  matrix  diagonalization  form  of 
symmetric  convolution  presented  here.  Thus  both  matrices  underlying  the  trigonometric  image 
reconstruction  filter  are  either  diagonal  or  approximately  diagonal  and  a  scalar  filter  results. 
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To  help  illustrate  how  a  filter  can  recover  an  estimate  of  an  object  from  its  distorted,  noisy 
version,  consider  the  imaging  scenario  in  Figure  1.  In  the  figure,  0(«,,«2)  represents  the  original 


object.  The  sequence  /i(n,,«2)  represents  the  point  spread  function  ofthe  system  blurring  the 
object.  The  scenario  adds  noise  represented  by  the  sequence  w(«i,«2)  to  produce  the  received 
data  sequence  The  recovery  filter  in  the  scenario  has  a 

two-dimensional  impulse  response,  /(ni,«2),  so  that  the  estimate  of  the  object  is  6{n^,n2)  = 

(«, ,  Wj ).  The  goal  of  image  reconstruction  is  to  find  the  impulse  response  of  the  re¬ 
covery  filter,  /(«i,«2).  which  produces  the  best  possible  estimate  0(n,,n2)- 

If  the  noise,  w(«i,«2),  is  negligible,  then  the  recovery  system  can  use  an  inverse  filter  to 
produce  an  estimate  ^(«,,«2)  =  ^(«i,n2)  which  recovers  the  object  exactly.  Otherwise  an  opti¬ 
mum  choice  for  a  recoveiy  filter  is  the  Wiener  filter.  Inverse  and  Wiener  filters  are  two  tradi¬ 
tional  types  of  linear  image  reconstruction  filters.  They  both  have  scalar  representations  in  the 
Fourier  domain.  Inverse  filters  simply  invert  the  Fourier  transform  coefficients  ofthe  degrada¬ 
tion  filter  to  recover  the  original  object  [15],  [27].  Wiener  filters  are  regularized  versions  of  in¬ 
verse  filters  designed  to  operate  in  the  presence  of  noise.  They  minimize  the  mean-squared  error 
between  the  estimates  they  recover  and  the  original  object  being  imaged  [21],  [38],  [44],  [52]. 
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This  research  shows  that  both  inverse  and  Wiener  filters  can  benefit  from  the  symmetric  convo¬ 
lution-multiplication  property  of  discrete  trigonometric  transforms.  The  results  derived  here 
create  new  versions  of  inverse  and  scalar  Wiener  filters  which  exist  in  the  transform  domain  of 
the  trigonometric  transforms. 

All  scalar  Wiener  filters,  including  the  traditional  Fourier  domain  and  new  trigonometric 
transform  domain  versions,  are  approximations  of  more  optimum  vector  or  generalized  Wiener 
filters  [38].  A  vector  Wiener  filter  accounts  for  the  off-diagonal  elements  present  in  the  trans¬ 
form  domain  representation  of  the  object  correlation  matrix.  Vector  Wiener  filters  provide  better 
estimates  than  scalar  Wiener  filters  [13],  [14],  but  they  are  more  computationally  intense  than 
scalar  Wiener  filters.  Pratt  [38]  shows  that  the  mean-squared  error  performance  of  vector  Wiener 
filters  in  the  transform  domain  is  independent  of  the  choice  of  transform.  Different  scalar  ap¬ 
proximations  of  a  vector  Wiener  filter  will,  however,  yield  different  levels  of  performance  based 
on  how  well  the  transform  domain  matrices  underlying  the  filter  approximate  diagonal  matrices. 
The  trigonometric  scalar  Wiener  filter  derived  here  yields  better  performance  than  its  Fourier 
counterpart  because  the  discrete  cosine  transform  yields  a  more  accurate  diagonal  approximation 
of  an  object  correlation  matrix  for  a  highly-correlated  image  than  a  discrete  Fourier  trans¬ 
form  [27]. 

Traditional  Fourier-domain  scalar  Wiener  filters  are  computationally  less  intense  than 
nonlinear  iterative  techniques,  but  suffer  from  poorer  mean-squared  error  performance.  The  fil¬ 
ters  which  result  from  this  research  require  fewer  computations  than  Fourier  scalar  Wiener  filters 
while  demonstrating  better  mean-squared  error  performance.  They  thus  provide  one  possible 
solution  to  the  Air  Force's  overall  problem  of  finding  image  reconstruction  techniques  which  are 
computationally  more  efficient  than  nonlinear  techniques,  but  which  have  better  performance 
than  existing  linear  methods. 


10 


This  research  effort  represents  the  first  time  the  symmetric  convolution-multiplication 
property  of  discrete  trigonometric  transforms  [34]  has  ever  been  applied  to  image  reconstruction. 
The  focus  is  therefore  on  applying  the  new  technique  to  the  most  basic  image  reconstruction 
methods.  To  apply  this  new  technique  to  create  new  versions  of  inverse  and  Wiener  filters  re¬ 
quires  a  few  basic  assumptions  which  must  underlie  the  model  of  the  basic  imaging  scenario. 
These  assumptions  involve  the  form  of  the  filters  which  cause  the  degradation,  the  statistics  of 
the  object  being  imaged,  and  the  noise  in  the  model. 

The  image  reconstruction  techniques  of  inverse  and  Wiener  filtering  are  linear  and  shift- 
invariant.  The  techniques  must  therefore  assume  that  the  degradation  is  caused  by  a  system 
which  is  also  linear  and  shift-invariant.  The  implication  of  linearity  and  shift-invariance  is  that 
the  form  of  the  point  spread  function  degrading  the  object  is  the  same  for  each  pixel  of  the  ob¬ 
ject.  This  assumption  is  valid  for  real-world  imaging  situations  [16].  The  trigonometric  inverse 
and  Wiener  filters  derived  here  require  knowledge  of  the  point  spread  function  of  the  degrading 
system.  Fourier  domain  inverse  and  scalar  Wiener  filters  also  require  a  known  point  spread 
function.  In  practice  the  point  spread  function  of  the  degrading  system  is  often  not  known  ex¬ 
actly.  In  many  cases  it  can  be  modeled  as  a  random  process  with  reasonable  assumptions  made 
for  its  statistics.  For  example,  if  the  distortion  in  an  image  is  caused  by  atmospheric  turbulence, 
then  reasonable  approximations  to  a  wide  variety  of  atmospheric  point  spread  functions  can  often 
be  found  [40]. 

An  additional  limitation  within  the  imaging  model  for  this  application  of  the  symmetric 
convolution-multiplication  property  is  related  to  the  statistics  of  the  object  being  imaged.  The 
object  statistics  must  be  assumed  to  be  wide-sense  stationary  and  highly  correlated.  A  wide- 
sense  stationary  object  has  a  constant  mean  and  an  autocorrelation  which  depends  only  on  the 
difference  between  the  locations  of  the  samples  in  an  image  [7].  If  the  vector  6  represents  the 
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object  sequence  e{n^,nj),  the  mean  of  the  vector  6  is  the  vector  fig.  For  a  wide-sense  station- 
aiy  object,  the  vector  Hg  must  first  be  constant  with  the  same  value  for  each  of  its  elements. 

The  second  condition  for  the  vector  9  to  be  wide-sense  stationaiy  is  that  the  w-nth  element  of  its 
correlation  matrix,  ,  must  only  depend  on  the  difference  between  m  and  n.  An  image 

with  highly-correlated  pixels  will  have  only  minor  fluctuations  in  the  values  of  its  pixels  in  a 
small  region.  Both  assumptions  of  wide-sense  stationarity  and  highly-correlated  pixels  are  nec¬ 
essary  for  the  correlation  matrix  to  become  almost  diagonal  in  the  trigonometric  transform  do¬ 
main.  A  good  scalar  approximation  of  the  filter  will  result  if  the  correlation  matrix  becomes  al¬ 
most  diagonal  in  the  transform  domain. 

The  imaging  model  also  includes  assumptions  regarding  the  noise  in  the  system.  The 
model  requires  that  the  noise  samples  be  zero-mean,  additive,  of  uniform  variance,  uncorrelated 
with  each  other,  and  independent  of  the  object.  These  assumptions  imply  that  the  amount  of 
noise  present  in  each  pixel  of  the  image  is  completely  unrelated  to  both  the  object  intensity  and 
the  amount  of  noise  present  in  any  other  pixel.  In  most  imaging  situations,  this  noise  model  is 
not  as  accurate  as  one  which  assumes  that  the  noise  is  an  object-dependent  Poisson  random  proc¬ 
ess  [17].  The  filtering  techniques  derived  here  do  not  use  a  more  accurate  Poisson  noise  model 
because  this  research  provides  a  first  look  at  the  benefits  of  symmetric  convolution  to  image  re¬ 
construction.  Earlier  Fourier  domain  inverse  and  Wiener  filters  use  exactly  the  same  noise 
model  as  that  assumed  here.  A  logical  next  step  for  future  research  would  be  to  extend  the  the¬ 
ory  developed  here  to  incorporate  object-dependent  photon  noise. 

The  assumptions  made  for  the  trigonometric  transform  domain  filters  derived  here  are  the 
very  same  assumptions  underlying  traditional  Fourier-domain  filters.  These  assumptions  allow 
the  performance  of  the  inverse  and  scalar  Wiener  filters  for  trigonometric  transforms  to  be  real- 
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istically  compared  to  their  traditional  Fourier  domain  ancestors.  The  research  performed  for  this 
dissertation  represents  the  first  time  the  symmetric  convolution-multiplication  property  of  trigo¬ 
nometric  transforms  has  ever  been  applied  to  image  reconstruction  of  any  type,  so  the  compari¬ 
son  to  earlier  techniques  is  warranted.  The  performance  results  of  these  brand  new  trigonometric 
transform  domain  filters  show  that  this  technique  appears  to  be  promising. 

This  research  achieves  its  overall  goal  of  applying  the  recently-developed  symmetric  con¬ 
volution-multiplication  property  of  the  discrete  trigonometric  transforms  [34]  to  the  traditional 
image  reconstruction  problems  of  inverse  and  Wiener  filtering.  The  results  extend  the  existing 
theory  behind  the  symmetric  convolution-multiplication  property  by  recasting  the  problem  into 
vector-matrix  form  and  then  demonstrating  how  discrete  trigonometric  transform  matrices  di¬ 
agonalize  matrices  which  represent  symmetric  convolution.  The  development  then  presents 
vector-matrix  forms  of  the  symmetric  convolution  property  for  multidimensional  asymmetric  se¬ 
quences  which  represent  the  most  general  type  of  sequences  encountered  in  image  reconstruc¬ 
tion.  The  filtering  of  multidimensional  asymmetric  sequences  is  then  possible  because  symmet¬ 
ric  convolution  is  equivalent  to  multiplication  in  the  transform  domain  for  each  of  the  underlying 
types  of  symmetry  in  an  asymmetric  image. 

The  research  presented  here  uses  the  newly-derived  vector-matrix  form  of  symmetric  con¬ 
volution  to  calculate  for  the  first  time  inverse  and  scalar  Wiener  filters  in  the  trigonometric  trans¬ 
form  domain.  The  new  forms  of  the  inverse  and  scalar  Wiener  filters  closely  resemble  their 
traditional  Fourier  domain  counterparts.  These  results  demonstrate  how  a  scalar  Wiener  filter 
provides  better  mean-squared  error  performance  for  symmetric  point  spread  functions  while  re¬ 
ducing  the  required  number  of  computations.  The  trigonometric  transform  domain  realizations 
use  fewer  calculations  because  the  trigonometric  transform  of  a  real  sequence  is  also  all  real. 

The  Fourier  transform  produces  complex  sequences  in  the  transform  domain. 
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The  remainder  of  this  dissertation  is  organized  as  follows.  Chapter  II  provides  some  back¬ 
ground  on  the  discrete  Fourier  transform  and  its  diagonalizing  forms,  inverse  and  Wiener  filter¬ 
ing,  and  the  symmetric  convolution-multiplication  property  of  discrete  trigonometric  transforms. 
The  notion  of  how  one  and  two-dimensional  discrete  Fourier  transform  matrices  diagonalize  cir- 
culant  and  block  circulant  matrices,  respectively,  serves  as  a  comparison  to  similar  diagonalizing 
results  derived  for  the  symmetric  convolution-multiplication  property  of  trigonometric  trans¬ 
forms.  The  background  on  inverse  and  Wiener  filtering  in  the  Fourier  domain  will  later  be  com¬ 
pared  to  new  results  for  these  filters  in  the  trigonometric  transform  domain.  The  background  on 
the  symmetric  convolution-multiplication  property  of  trigonometric  transforms  is  essential  to 
understanding  how  this  property  is  later  extended  theoretically  and  then  applied  to  the  inverse 
and  Wiener  filtering  problems. 

All  of  the  material  presented  in  Chapter  II  summarizes  previously  conducted  work,  while 
the  material  in  Chapters  III  -  V  represents  the  results  of  new  work  performed  for  this  dissertation. 
In  Chapter  HI,  the  symmetric  convolution-multiplication  property  of  discrete  trigonometric  trans¬ 
forms  is  derived  using  vector-matrix  methods  and  then  extended  to  multidimensional  and  asym¬ 
metric  sequences.  Chapter  IV  presents  the  derivation  of  one  and  two-dimensional  inverse  filters 
expressed  in  the  trigonometric  transform  domain  of  symmetrically  convolved  sequences.  Also  in 
Chapter  IV  are  derivations  of  one  and  two-dimensional  trigonometric  transform  versions  of  sca¬ 
lar  Wiener  filters.  Chapter  V  provides  an  example  of  recovering  a  distorted  image  using  inverse 
and  scalar  Wiener  filters  and  details  the  mean-squared  error  performance  of  the  scalar  Wiener 
filter.  A  conclusion  appears  as  Chapter  VI. 
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II.  Background 


This  chapter  provides  background  on  the  discrete  Fourier  transform  and  its  diagonalizing 
forms,  inverse  and  Wiener  filtering,  and  the  symmetric  convolution-multiplication  property  of 
discrete  trigonometric  transforms.  The  concept  of  how  one  and  two-dimensional  discrete  Fourier 
transform  (DFT)  matrices  diagonalize  circulant  and  block  circulant  matrices  will  later  serve  as  a 
comparison  to  similar  diagonalizing  results  for  the  symmetric  convolution-multiplication  prop¬ 
erty  of  trigonometric  transforms.  In  the  second  section,  inverse  and  Wiener  filters  are  derived  in 
the  Fourier  domain.  The  traditional  representation  of  these  filters  will  be  similar  to  new  trigo¬ 
nometric  transform  versions  derived  in  later  chapters.  The  final  section  of  this  chapter  intro¬ 
duces  the  symmetric  convolution-multiplication  property  of  trigonometric  transforms.  This 
chapter  provides  a  baseline  from  which  the  property  is  later  extended  to  a  more  general  class  of 
signals  and  then  used  to  derive  inverse  and  Wiener  filters  in  the  trigonometric  transform  domain. 

2.1  Diagonalizing  Forms  of  the  Discrete  Fourier  Transform 

The  well-known  circular  convolution-multiplication  property  of  the  DFT  [37]  for  one  and 
two  dimensions  is  reviewed  in  this  section.  One  and  two-dimensional  DFT  matrices  are  defined 
and  then  shown  to  diagonalize  circulant  and  block  circulant  matrices,  respectively. 

Hunt  [23]  was  the  first  to  observe  that  a  DFT  matrix  diagonalizes  a  circulant  matrix  which 
performs  circular  convolution.  Consider  two  finite  one-dimensional  sequences,  h{n)  and  6{ri), 
which  equal  zero  outside  the  interval  0  <  «  <  iV  - 1.  The  circular  convolution  of  the  two  se¬ 
quences  is  expressible  in  vector-matrix  notation  as  d  =  h©0-  H^^d,  where  ©  represents  circu¬ 
lar  convolution.  The  NxN  matrix  is  circulant  as  indicated  by  the  subscript  'C,  and  has  the 
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vector  h  as  its  first  column.  Hunt  proves  that  is  transformed  into  a  diagonal  matrix  through 
the  relation 

Hc  =  W,^,W^\  (1) 

In  Eq.  (1),  is  the  NxN  DFT  matrix  with  m-nth  entry 

[^w'L„  m,n=  0,  1,  (2) 

Its  inverse,  ,  is  the  NxN  inverse  DFT  matrix  with  m-nth  entry 

Wn\  =— exp{+7^wn};  m,n=  0,  I,  ...,N-\.  (3) 

The  NxN  matrix  in  Eq.  (1)  is  diagonal  with  the  elements  of  the  vector  along  the 
diagonal.  Hunt  [23]  shows  that  the  elements  of  are  the  eigenvalues  of  The  boldface 
script  notation  indicates  a  matrix  in  the  transform  domain,  while  the  subscript  'F'  indicates  the 
Fourier  transform  domain.  Thus  =  ^pW^^d  or  equivalently 

(4) 

where  ©  represents  a  point-wise  or  Hadamard  product  [22].  Equation  (4)  is  the  one-dimensional 
circular  convolution-multiplication  property  for  the  DFT  expressed  in  vector-matrix  form. 

Hunt's  results  extend  naturally  to  two  dimensions  [24]  using  Kronecker  (or  tensor)  prod¬ 
ucts  [18].  The  Kronecker  product,  0,  of  an  Mj  x  matrix,  A,  and  an  M2  x  N2  matrix,  B,  is 

14o«  l^lw-i)-» 

where  [^]„„  is  the  m-nth  entry  of  the  matrix  A.  The  Kronecker  product  A0B  has  dimension 
M,M2  ^  ^1^2-  I**  dimensions,  the  lexicographically-ordered  vector  0  =  vec{0}  provides 
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an  alternate  means  to  express  an  N^  x  N2  matrix,  0.  The  vec{«}  operation  converts  the  matrix 
into  an  N^N2xl  column  vector  of  stacked  transposed  rows  [24].  The  vector-matrix  expression 
for  the  two-dimensional  DFT  is  then 

=  (6) 

Note  that  it  requires  fewer  calculations  to  compute  0p  =  where  <0p  is  the  matrix 

representation  of  the  lexicographically-ordered  vector,  &p.  This  matrix  form  exists  because  the 
Fourier  transform  acts  on  the  rows  and  columns  of  the  matrix  0  separately  [27].  The  connection 
between  the  one  and  two-dimensional  DFT  is,  however,  easier  to  visualize  from  Eq.  (6).  Equa¬ 
tion  (6)  is  also  the  only  way  to  calculate  the  second  order  moments  of  the  vector  6  in  the  trans¬ 
form  domain  when  9  is  later  modeled  as  a  random  vector. 

To  extend  his  one-dimensional  results  to  two  dimensions.  Hunt  [24]  defines  the  discrete 
two-dimensional  circular  convolution  of  two  iVj  x  N2  matrices,  H  and  0,  as  a  vector-matrix 
operation.  The  matrices  FT  and  0  represent  finite  two-dimensional  sequences  /?(«,, «2)  and 
0(«j,«2)  which  are  zero  outside  the  region  0  <  - 1,  0  <  Wj  ^  ^2  “  1-  The  vector-matrix 

form  of  two-dimensional  circular  convolution  is  </  =  Hp(-0,  where  d  and  9  are  lexicographi¬ 
cally-ordered  vectors.  The  matrix  Hgc  is  an  N^N2  x  block-circulant  matrix,  which  is  a 
matrix  of  partitioned  blocks  of  circulant  matrices  arranged  in  a  circulant  pattern.  The  first  col¬ 
umn  of  the  matrix  will  be  the  vector  h  =  Hunt  shows  that  the  N1N2  x  N1N2  two- 

dimensional  DFT  matrix,  diagonalizes  Hgc-  Specifically, 

Hbc  =  K  ®  ®  <).  (7) 
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where  is  diagonal  with  the  vector  0  along  the  diagonal.  Equation  (7) 

implies  that 

«  ®  V  =  K‘  ®  ©  (fv;;;  0  (8) 

which  is  the  two-dimensional  version  of  Eq.  (4). 

The  circular  convolution-multiplication  property  of  DFTs  greatly  reduces  the  number  of 
calculations  in  image  reconstruction  problems  which  involve  linear  filtering.  The  computational 
cost  savings  in  implementing  Eq.  (8)  instead  of  d  =  arises  because  it  requires  fewer  com¬ 
putations  to  compute  the  transforms,  point-multiply  the  results,  and  then  apply  an  inverse  trans¬ 
form  than  it  does  to  compute  the  large  vector-matrix  multiplication  directly.  Note  that  the  DFT  is 
not  a  unique  transform  for  implementing  filtering  operations  in  the  transform  domain.  This  re¬ 
view  of  the  diagonalizing  forms  for  the  DFT  matrix  has  established  a  reference  point  from  which 
to  analyze  and  later  apply  the  results  of  more  recent  developments  in  the  diagonalizing  forms  of 
the  trigonometric  transforms. 

2.2  Signal  and  Image  Reconstruction  in  the  Fourier  Domain 

The  vector-matrix  diagonalizing  forms  of  the  previous  section  play  an  important  role  in 
deriving  traditional  Fourier  domain  inverse  and  Wiener  filters.  These  filters  can  reconstruct  one¬ 
dimensional  signals  and  two-dimensional  images  which  are  degraded  by  linear  processes.  This 
background  material  on  inverse  and  Wiener  filtering  in  the  Fourier  domain  will  later  serve  as  a 
comparison  to  similar,  but  improved  versions  of  these  filters  derived  in  the  trigonometric  trans¬ 
form  domain. 

2.2.1  Inverse  filtering.  The  derivations  of  inverse  filters  for  one  and  two  dimensions  in 
the  Fourier  domain  are  presented  in  this  subsection.  The  one-dimensional  circular  convolution- 
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multiplication  property  in  Eq.  (4)  is  valid  for  any  x  1  vectors  d,  h,  and  0,  which  represent  N- 
point  sequences  d(n),  h(n),  and  0{n).  All  three  sequences  equal  zero  outside  the  interval 
0<n<  N -1.  In  one-dimensional  signal  reconstruction,  the  vectors  and  sequences  take  on  spe¬ 
cial  meaning  beyond  that  of  just  being  arbitrary  sequences  for  which  the  circular  convolution 
d  =  HqO  holds.  The  sequence  h{n)  represents  the  finite  impulse  response  of  a  linear  shift- 
invariant  filter  which  distorts  a  signal  d{ri)  to  produce  a  data  sequence  d(n).  The  goal  of  signal 

A. 

reconstruction  is  to  obtain  an  estimate,  ^(n),  of  9{n)  from  the  corrupted  data  sequence  d{n). 

A  notational  variation  of  this  problem  is  to  replace  the  above  sequences  with  vectors.  The 
goal  of  the  problem  is  to  find  the  impulse  response  of  a  filter,  represented  by  the  vector /,  which 
recovers  the  vector  0  from  d,  given  knowledge  of  h.  Because  the  underlying  form  of  convolu¬ 
tion  for  the  DFT  is  circular,  the  vector / can  be  represented  as  a  circulant  matrix,  F^,  so  that 

6  =  F(,d  =  FcH^d.  The  matrix  Fq  will  recover  the  vector  6  exactly  if  ,  or  equiva¬ 

lently  in  the  Fourier  domain  if 

(9) 

It  follows  that  Both  of  the  matrices  are  diagonal,  and  can  be  completely 

characterized  by  their  elements  &nd‘t¥p{k),  so  that  Eq.  (9)  reduces  to  the  scalar  equation 


— 


1 

l^pik)  ’ 


(10) 


for  ik  =  0,  1 ,  . . . ,  - 1,  provided  ‘!^p{k)^Q.  Recall  that  script  letters  denote  transform  domain 

quantities,  and  the  subscript  'F'  denotes  the  Fourier  domain.  Equation  (10)  is  the  one¬ 
dimensional  inverse  filter  expressed  in  the  Fourier  domain. 


19 


The  two-dimensional  imaging  scenario  uses  the  NiN2  x  1  vectors  d  and  0  to  represent  a 
detected  image  and  an  original  object,  respectively.  The  vectors  are  lexicographic  reorderings  of 
the  X  N2  matrices  D  and  0.  The  two-dimensional  finite  impulse  response,  H,  of  a  system 
which  distorts  an  object  is  its  point  spread  function  (PSF).  The  lexicographically-ordered  vector 
h  =  vec{ir}  represents  the  PSF.  As  in  the  previous  section,  the  vector  h  is  the  first  column  of 
the  N^N2  X  AjA/j  block  circulant  matrix  In  the  two-dimensional  circular  convolution  re¬ 
lation,  d  =  the  matrix  blurs  the  object  to  produce  the  distorted  data. 

Just  as  in  one  dimension,  this  classical  image  reconstruction  problem  is  to  find  the  two- 
dimensional  impulse  response  of  a  filter  which  recovers  an  estimate  6  from  d,  given  h.  The  im¬ 
pulse  response  is  represented  by  the  lexicographically-ordered  vector /.  The  two-dimensional 
convolution  operation  of  the  filter  is  implemented  as  a  block  circulant  matrix,  F^f.,  with  the 

vector / as  its  first  column.  The  expression  to  recover  the  estimate  then  becomes  9 

Fg^Hg^d-  The  object  vector  6  is  exactly  recovered  if  ^BC  =  Hg^.  This  problem  is  commonly 
presented  in  the  Fourier  domain  [27],  [15]  as 

K  ®  wv.  K‘  ®  f'ii) = [K  ®  ®  w”;;.')]"'-  (11) 

From  Eq.  (1 1)  it  follows  that  “Pp  =  where  both  matrices  are  diagonal  and  can  be  represented 
by  their  diagonal  elements  pp  (A:, ,  )  and  ( A, ,  A2 ).  The  scalar  form  of  Eq.  ( 1 1 )  is  thus 


Pp(Ji\,k2^ 


1 

P^p(k^,k2) 


(12) 


for  A,=0,  1,  ...,^"1-1;  A2=0,  1,  ...,N2-l;  provided  ‘»p(k^,k2)^0.  Equation  (12)  is  the  two- 
dimensional  Fourier  domain  inverse  filter.  Note  its  similarity  to  Eq.  (10). 
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One  limitation  of  inverse  filtering  in  the  Fourier  domain  [27],  [15]  is  that  for  most  praeti- 
cal  realizations,  takes  on  very  small  values  for  increasing  k^  and  k2.  These  small 

values  produce  very  high  gains  in  7p{k^,k2)  at  the  higher  spatial  frequencies  in  Eq.  (12).  This 
problem  becomes  quite  serious  when  the  image  model  expands  to  incorporate  noise  with  a  uni¬ 
form  power  spectral  density  across  all  frequencies.  The  solution  to  this  problem  involves  regu¬ 
larizing  the  high  frequency  gain  of  the  inverse  filter. 

2.2.2  Wiener  filtering.  The  development  of  scalar  Wiener  filters  for  one  and  two  dimen¬ 
sions  in  the  Fourier  transform  domain  is  reviewed  in  this  subsection.  Wiener  filters  bear  the 
name  of  Norbert  Wiener  who  pioneered  them  in  the  1940s  during  studies  on  time  series  analy¬ 
sis  [52].  Helstrom  [21]  and  Slepian  [44]  later  applied  Wiener  filters  to  the  image  reconstruction 
problem.  These  filters  introduce  a  degree  of  regularization  to  the  inverse  problem  by  minimizing 
the  mean-squared  error  between  the  object  estimate  they  produce  and  the  original  object.  They 
have  better  mean-squared  error  performance  than  inverse  filters  when  reconstructing  noisy  sig¬ 
nals  or  images. 

The  data  model  of  the  previous  sections  must  have  an  added  term  to  incorporate  noise  so 
that  the  model  now  becomes  d  =  H(.0  +  w  for  the  one-dimensional  circulant  case.  In  the  new 
model,  iv  is  an  N  x  \  zero-mean  uniform-variance  noise  vector  whose  samples  are  uncorrelated 
both  with  the  object  vector,  6,  and  with  each  other.  The  model  assumes  that  the  object  vector, 

6,  is  also  a  random  vector  which  will,  in  general,  have  a  mean  vector  Hq  =  E[d^  =  9.  The  op¬ 
erator  £{•}  and  the  overbar,  •  ,  are  equivalent  ways  to  denote  expectation.  The  vector  9  will 
also  have  a  covariance  matrix  ^ee  -  E^9-  Under  these  assumptions,  Wiener 

filters  produce  an  estimate  represented  by  the  vector  9  which  minimizes  the  expected  value  of 
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the  square  of  the  2-norm  of  the  error  vector,  e  =  6-9.  If  denotes  the  2-nonn  of  the  vector 

s,  then  the  linear  minimum  mean-squared  error  estimate,  occurs  at  minjllfll^l  or  equiva¬ 

lently  minjf^f  . 

Solutions  to  the  problem  of  finding  filters  which  recover  9  from  d  under  these  conditions 
appear  in  many  statistical  signal  processing  texts  [28],  [43],  [48].  The  solution  which  follows  is 
from  [28]  and  is  the  linear  minimum  mean-squared  error  estimate  which  results  from  the  Baye¬ 
sian  Gauss-Markov  Theorem  under  the  conditions  described  above: 

9^fie  +  CeeHl[HcC^l  +  cJi\d-HcMe].  (13) 

If  the  mean  of  the  vector  9  is  constant,  then  without  loss  of  generality  [28],  the  mean  can  be  as¬ 
sumed  to  be  the  zero  vector,  0,  so  that  Eq.  (13)  becomes 

9  =  ReM[HcReeHl  +  d.  (14) 

In  Eq.  (14)  the  covariance  matrices  of  Eq.  (13)  are  replaced  with  correlation  matrices  defined  by 
Rgg  =  00^ .  In  this  case,  Cgg  will  equal  Rgg  because  of  the  assumption  that  the  vector  Hg  -  0. 

It  is  not  difficult  to  derive  Eq.  (14)  using  standard  vector-matrix  minimization  techniques 
as  outlined  in  [38].  First  observe  that  the  goal  of  the  problem  is  to  find  the  filter  which  mini¬ 
mizes  ll^llj  =  s  =Tr|f^£'|  =  Trjfi’f^  .  The  operation  Tr{«}  denotes  the  trace  of  a  matrix. 

The  transpose  denoted  by  the  superscript  T'  is  used  because  it  is  assumed  that  all  sequence  do¬ 
main  quantities  are  real.  Substituting  9-9  for  the  error  vector,  e,  produces 

ii4=Tr{£{(«-9)(«-»)''}} 

=  Tr^E\99'^  -99'^  -99'^  +  99'^'^  .  .  . 
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=  Tr  j  £  {6>(9^  -  2(96»^  +  Jj . 


(15) 


The  intermediate  steps  of  Eq.  (15)  use  the  fact  that  for  any  two  real  vectors,  Tr|flA^|  = 

Tr|6^a|  =  b^a  =  a^b  =  Tr|a^A|  =  TrjAa^j.  If  is  the  circulant  matrix  version  of  the  one¬ 
dimensional  impulse  response  of  the  reconstruction  filter,  then  substituting  for  6  =  F(.d  = 

F(.Hcd  +  F(.w  yields 

114  - Tri^E^ee^ -2{FcHcd  +  Fcw)0^  +(Fc^c^  +  f^c»^)(^c^c^+^c»^)^)}-  (16) 

After  expanding  and  recognizing  that  because  the  object  and  noise  are  uncorrelated,  wd^  = 

6w^  =  0,  an  A  X  iV  zero  matrix,  Eq.  (16)  becomes 

\\el=Tr[Ree-2FcHcRee  +  Fc[HcReeHl^RjjF^],  (17) 

where  Rgg  =  69^  and  =  ww^ .  The  next  step  in  finding  the  filter,  F^,  which  minimizes 
Ifl^,  is  to  take  the  first  derivative  of  Eq.  (17)  with  respect  to  which  yields 

^-^  =  -2R,^l  +  2fJHcR,oH^c  +^w)-  (18) 

^Fc  ^  ’ 

Taking  the  derivative  of  Eq.  (17)  requires  using  the  matrix  derivatives  <7Tr|x4  |  jdX  =  and 

dTx^XAX^^jdX  ^x{A->rA^^  =2X4,  if  A  is  symmetric  [54]. 

Setting  Eq.  (18)  equal  to  the  zero  matrix  and  solving  for  F^  produces  the  filter  which 
minimizes  the  mean-squared  error  between  the  original  object  and  its  estimate, 

Fc  =  ReeHl[HcReeHl^Kj[\  (19) 
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Equation  (19)  is  exactly  equivalent  to  Eq.  (14)  because  0  =  F(^d.  Converting  Eq.  (19)  to  the 
Fourier  domain  using  Eq.  (1),  produces 


-1-1 


(20) 


The  superscript  'H'  which  denotes  the  Hermitian  or  conjugate  transpose  is  now  required  because 
Fourier  domain  quantities  are  complex.  Solving  Eq.  (20)  for  7^  and  distributing  the  Hermitian 
operator  results  in 


where  1?^  =  because  is  diagonal.  Equation  (21)  reduces  to 


(21) 


^F  -  R0^0f,'^F^FReF0p‘^F  RltiFKi^\  ’  (22) 

where  ==  Ree^N  RwfWf  ~  RwvWn  •  Equation  (22)  is  the  general  or  vector 

Wiener  filter  in  the  Fourier  transform  domain  [38]. 

Note  that  in  Eq.  (22)  the  matrices  Tr,  and  are  diagonal  by  definition,  and  the 
matrix  is  diagonal  because  of  the  assumption  of  uniform-variance  uncorrelated  noise 
samples.  The  matrix  R@^@^  will  be  approximately  diagonal  under  the  assumption  of  wide-sense 

stationarity  for  the  object.  The  approximation  arises  because  the  discrete  Fourier  transform  does 
not  exactly  diagonalize  the  symmetric  Toeplitz  form  of  the  correlation  matrix  for  a  wide-sense 
stationaiy  object  [48].  The  scalar  equation  which  results  from  approximating  Rq^^  by  retain¬ 
ing  only  its  diagonal  elements  is  then 
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The  terms  TOpik)  and  0\{k)  are  the  diagonal  elements  of  and  and  represent  the 

power  spectral  densities  of  the  noise  and  object  respectively.  The  terms  and 
the  elements  of  the  diagonal  matrices  and  ^p.  Equation  (23)  also  uses  the  fact  that 


‘!i‘p{k)'i^*p{k)  =  The  Fourier  domain  filter  in  Eq.  (23)  is  referred  to  as  the  one¬ 

dimensional  scalar  Wiener  filter.  Under  the  assumption  of  wide-sense  stationarity  for  the  object, 
it  is  a  good  approximation  to  the  vector  Wiener  filter  [38]. 

The  extension  of  the  preceding  development  to  two  dimensions  is  straightforward.  The 
data  model  incorporating  noise  becomes  d  =  Hp^G+w.  The  N^N2  x  #,^2  block  circulant  deg¬ 
radation  matrix  has  the  lexicographically-ordered  vector  ft  -  vec{JT}  which  represents  the 
PSF  as  its  first  column.  The  vectors  d  and  6  are  N.^N2'x\  lexicographically-ordered  vectors 
representing  the  data  and  the  original  object  matrices  D  and  0  respectively.  The  lexicographi¬ 
cally-ordered  vector  w  is  an  N^N2  x  1  zero-mean  uniform-variance  noise  vector  whose  samples 
are  again  assumed  uncorrelated  both  with  the  object  vector,  G,  and  with  each  other.  The  object 
vector,  6,  is  again  random.  It  has  a  constant  mean  vector  which  equals  0  without  loss  of 


generality.  The  autocorrelation  matrix  of  the  vector  6  is  ^ee  =  66^.  The  linear  minimum 


mean-squared  error  estimate  represented  by  the  vector  0  is  then  [28] 


-1 


d. 


(24) 


Q  —  Ree^ BC^ee^ Bc 
The  only  differences  between  Eqs.  (14)  and  (24)  are  that  in  the  two-dimensional  case  the  degra¬ 
dation  matrix,  Hgc,  is  block  circulant,  the  autocorrelation  matrix,  Rgg,  is  block  symmetric 


Toeplitz,  and  the  vectors  0,  6,  h,  w,  and  d  are  all  lexicographically-ordered.  From  Eq.  (24)  it 
follows  that  the  two-dimensional  impulse  response  of  the  reconstruction  filter  is 


^Bc  ~  ^ee^BC 


^ BC^ee^ BC  •  (25) 

Expanding  Eq.  (25)  using  Eq.  (7)  produces 

7f  ~  ^epBp^F  ^  ’  (26) 

in  the  Fourier  domain.  As  in  the  one-dimensional  case,  the  vector  Wiener  filter  of  Eq.  (26)  is 
well-approximated  by  the  scalar  equation  [15] 


^f(^1>^2) 


6>f(^1.*2) 


(27) 


which  retains  only  the  diagonal  elements  of  the  matrices  in  Eq.  (26).  The  terms  7(/p{k^,k2)  and 


0^(^l,^2)  are  the  diagonal  elements  of  and  ^0p@p  which  represent  the  two-dimensional 

power  spectral  densities  of  the  noise  and  object  respectively.  The  Fourier  domain  filter  in 
Eq.  (27)  is  the  two-dimensional  scalar  Wiener  filter  for  reconstructing  a  corrupted  image  in  noise 
given  knowledge  of  the  degrading  PSF  which  distorted  it. 

The  overall  goal  of  this  research  is  to  derive  new  expressions  for  the  one  and  two- 
dimensional  inverse  and  scalar  Wiener  filters  of  Eqs.  (10),  (12),  (23),  and  (27)  which  lie  not  in 
the  Fourier  domain,  but  which  instead  lie  in  the  trigonometric  transform  domain.  These  filters 
will  also  employ  the  newly-developed  symmetric  convolution-multiplication  property  of  trigo¬ 
nometric  transforms  [34].  To  derive  expressions  for  these  filters  in  the  trigonometric  transform 
domain  which  employ  this  new  property,  it  is  necessary  to  first  provide  some  background  on  the 
symmetric  convolution-multiplication  property  itself. 
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2.3  Symmetric  Convolution  and  the  Discrete  Trigonometric  Transforms 

This  section  provides  background  on  the  recently-developed  symmetric  convolution- 
multiplication  property  of  discrete  trigonometric  transforms  (DTTs)  [34].  The  material  presented 
here  serves  as  an  essential  baseline  from  which  the  property  will  be  theoretically  extended  to  a 
broader  class  of  signals  in  later  chapters  and  then  applied  to  traditional  inverse  and  Wiener  filter¬ 
ing  problems. 

The  discrete  cosine  transform  was  first  introduced  in  1974  [1].  Since  then  it  has  been  ex¬ 
panded  into  an  entire  family  of  trigonometric  transforms  [26],  [39]  consisting  of  sixteen  DTTs 
which  are  even  and  odd-length  versions  of  the  discrete  sine  and  cosine  transforms  (DSTs  and 
DCTs).  The  family  lies  at  the  heart  of  many  image  transform  coding  applications  [27],  [33], 

[39].  Although  useful  for  transform  coding,  the  DTTs  have  not  proved  veiy  useful  in  a  wide  va¬ 
riety  of  image  filtering  applications,  because  until  recently  no  convolution-multiplication  prop¬ 
erty  existed  for  the  entire  family. 

Martucci  [34]  has  recently  developed  a  convolution-multiplication  property  for  the  family 
of  DTTs.  He  defines  symmetric  convolution  as  the  form  of  convolution  for  DTTs.  His  results 
are  analogous  to  the  circular  convolution-multiplication  property  of  the  DFT  [37].  The  symmet¬ 
ric  convolution-multiplication  property  states  that  an  inverse  trigonometric  transform  of  the 
product  of  the  trigonometric  transforms  of  two  sequences  yields  the  same  result  as  the  symmetric 
convolution  of  the  two  sequences  [34].  The  symmetric  convolution-multiplication  property  exists 
for  forty  different  combinations  of  the  sixteen  transforms  in  the  DTT  family  based  on  the  under¬ 
lying  symmetric  periodicities  of  the  different  DTTs. 

There  are  four  ways  to  symmetrically  extend  a  finite  sequence  about  a  single  point  of 
symmetry.  These  are  whole-sample  symmetry  (WS),  whole-sample  antisymmetry  (WA),  half- 
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sample  symmetry  (HS),  and  half-sample  antisymmetry  (HA).  An  example  of  each  appears  in 
Figure  2.  Note  that  the  point  of  symmetry  for  the  WS  sequence  is  the  end  point  in  the  finite  se¬ 


quence  before  extension,  and  the  point  of  symmetry  for  the  WA  sequence  is  a  zero  which  must 
appear  after  the  end  point  before  extension.  The  points  of  symmetry  for  both  the  HS  and  HA 
sequences  lie  one-half  sample  beyond  the  end  points  in  each  finite  sequence  before  extension. 

There  are  16  symmetric  periodic  sequences  (SPS's)  which  result  from  symmetrically  ex¬ 
tending  a  finite  sequence  to  the  left  using  one  of  the  four  ways  and  to  the  right  using  possibly  a 
different  way.  A  convention  for  naming  each  of  the  16  SPS's  is  to  label  first  the  left  symmetric 
extension  and  then  the  right  symmetric  extension.  For  example  a  WSHA  sequence  would  exhibit 
whole-sample  symmetry  to  the  left  and  half-sample  antisymmetry  to  the  right. 

A  finite  sequence  is  converted  into  an  SPS  by  applying  a  symmetric  extension  operator. 

For  example,  the  WSHA  extension  of  x{n),  for«  =  0,  1,  ...,  N-\,  is 


x{ri)  =  e 


WSHA 


x(n), 

-xiM-n), 


n  =  0,  1,  ...,  N -\ 
n  =  N,  ...,  M-\, 


(28) 


where  M  =  2N-\.  The  symmetric  extension  operator  {•}  is  one  of  sixteen  different 
symmetric  extension  operators  ~  one  for  each  type  of  SPS.  If  jc  is  a  column  vector  whose  indi- 
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vidual  entries  are  x{n),  then  Eq.  (28)  is  equivalent  to  the  vector-matrix  equation,  x  = 


where  x  is  an  Mxl  column  vector.  The 

Mx  N  matrix  is 

defined  by 

X(0) 

’1 

x(0) 

x(l) 

1 

x(l) 

x(2) 

1 

jc(2) 

X  -  Ei^shaX  - 

jc(V-l) 

= 

1 

_x(N-l)_ 

-x(N-l) 

-1 

-x(2) 

-1 

-^(1) 

0  -1 

The  definitions  of  the  remaining  symmetric  extension  operators  follow  in  a  similar  fashion. 

Each  of  the  16  DTTs  has  an  underlying  symmetry  for  both  the  input  sequence  and  the  out¬ 
put  sequence  in  the  trigonometric  transform  domain.  This  behavior  is  similar  to  the  underlying 
periodicity  of  the  input  and  output  sequences  for  the  DFT.  The  underlying  input  and  output 
symmetries  for  each  DTT  as  well  as  its  form  and  other  information  appears  in  Table  1.  The  sub¬ 
scripts  on  each  transform  denote  the  type  (from  types  I  -  IV)  and  whether  the  transform  is  even  or 
odd.  Note  that  unlike  the  previous  notational  convention,  no  subscript,  'N,'  indicating  the  length 
of  the  transform  appears  in  the  table.  Omitting  the  subscript  preserves  space  in  the  table,  and  this 
information  appears  in  the  columns  for  the  input  and  output  ranges.  Also  the  constant  k„  which 
appears  in  the  definition  of  the  transforms  should  not  be  confused  with  the  one-dimensional 
transform  domain  index,  k,  or  with  the  two-dimensional  transform  domain  indices,  it,  or  itj. 

The  form  of  each  transform  in  Table  1  is  slightly  different  from  traditional  representations 
of  the  DTTs  [26],  [39],  [50].  Martucci  [34]  refers  to  the  form  in  Table  1  as  the  convolution  form 
of  the  DTTs  and  relates  each  transform  to  its  representation  in  terms  of  the  generalized  discrete 
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Table  1.  Forward  Discrete  Trigonometric  Transforms 


H 

m-n\h  Entry 

1  Input 

1  Output 

Equivalent 

Transform 

Symmetry 

Range 

Symmetry 

Range 

6u 

KL  =  2*.cc 

wsws 

o-^w 

wsws 

o-^w 

^O.O^WSfFS 

S,e 

KL=2sin( 

7imn\ 

If) 

WAWA 

1  -^N-l 

WAWA 

1  -^W-l 

J^0,0^fK4fVA 

KL  =  2cos| 

<  N  ) 

HSHS 

o->m 

WSWA 

0^  W-l 

^qI^HSHS 

m 

[‘52.L  =  2sin[ 

mi(n  + 1)^ 

N  ) 

HAHA 

0->N-l 

WAWS 

1->N 

[C3.L  =  2*,cc 

WSWA 

0->N-l 

HSHS 

0-^N-l 

^lO^WSfVA 

[‘^3eL  = 

WAWS 

1->N 

HAHA 

O^W-1 

j^lfi^WAWS 

KL  =  2«>s( 

) 

HSHA 

0->N-l 

HSHA 

0-^W-l 

01L^HSHA 

^4e 

[‘^4.L  =  2s'n( 

^(w  +  i)(«  +  i)) 

N  J 

HAHS 

0->N-l 

HAHS 

0-^W-l 

J^ll^HAHS 

(?lo 

HH 

WSHS 

O^N-1 

WSHS 

0->N-l 

0Ofi^WSHS 

[5,L  =  2sin(- 

2mnn 

2N-lJ 

WAHA 

WAHA 

1  ->W-1 

j^Qfi^WAHA 

&2o 

[c4„„  =  2/„co 

/2mn(»  +  i)'\ 

\  2N-1  J 

HSWS 

0-^N-l 

WSHA 

O^W-1 

^qI^HSWS 

B 

KL  =  2sin[ 

2mn(n  + 

2A^-1  J 

HAWA 

O-^N-2 

WAHS 

1  ->W-1 

J^qV^HAWA 

KL  =  2*„cc 

/2>r(m  +  i)n'\ 

V  2iy-i  J 

WSHA 

0-^N-l 

HSWS 

O^W-1 

[■^^L  =  2sjn(- 

2fr(m  -f  })n\ 

2N-1  J 

WAHS 

l->N-l 

HAWA 

O-^W-2 

j^lfi^WAHS 

<?4o 

KL  =  2cos( 

'2;r(m  +  i)(n  +  i)^ 

2N-1  J 

HSWA 

O^N-2 

HSWA 

O^W-2 

B 

KL  =  2/,sin 

f2^(m  +  i)(n  +  m 

1  2N-1  J 

HAWS 

0-^N-l 

HAWS 

O^W-1 

J011^HAWS 

ri/2,  n  =  0oTN  ri,  n  =  0,l,  ...,N-2 

"  [I,  n  =  l,  2,  N-l  "  [1/2,  n=N-l 


Fourier  transform  (GDFT)  [3].  The  GDFT  operator,  represented  by  the  script  character,  ^  {•}, 
is  expressible  as  a  matrix  with  m-«th  entry 

=exp(^-y^('«  +  «)(«+*)]i  m,n  =  Q,  1,  ...,  M-\.  (30) 
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The  boldface  character  in  Eq.  (30)  denotes  the  matrix  which  represents  the  operator, 

{•}.  The  subscripts  'a'  and  V  indicate  the  amounts  of  the  shifts  in  the  sequence  and  trans¬ 
form  domains.  The  subscript  'Af  denotes  the  Mx  M  dimension  of  the  transform  matrix  which 
results  from  extending  the  original  sequence  length  so  that  M  =  2N  for  N  even  or  M  =  2N-\ 
for  odd.  Expressed  in  this  form,  each  transform  in  Table  1  is  equivalent  to  the  GDFT  of  a 
symmetrically  extended  version  of  the  finite  input  sequence  with  the  same  symmetry  that  under¬ 
lies  the  trigonometric  transform.  For  example,  where  h(n)  is  a 

finite  sequence  which  the  type-I  DCT  assumes  has  WSWS  underlying  symmetry.  These  equiva¬ 
lent  relationships  also  appear  in  the  last  column  of  Table  1  for  each  transform.  The  values  of  a 
and  b  in  the  GDFT  for  these  cases  are  always  either  0  or 

Table  2  lists  similar  information  to  Table  1  for  each  inverse  DTT.  The  second  column  of 
Table  2  lists  the  relationship  between  the  forward  and  inverse  transforms. 

Based  on  the  underlying  principles  discussed  above,  Martucci  [34]  defines  the  symmetric 
convolution  of  two  sequences,  0(n)  and  h(n),  as 

^(«)  ==  .  (31) 

where  Si(n)  is  a  rectangular  window  of  length  L  which  extracts  the  representative  samples  from 
the  result.  The  operations  s^{6{n)]  and  S2{h{n)]  create  symmetric  periodic  extensions  of  the 

two  sequences  being  convolved.  The  convolution  operation  '©'  represents  either  circular  or 
skew-circular  convolution.  The  need  to  perform  skew-circular  convolution  arises  because  in  half 
of  the  forty  cases  of  symmetric  convolution,  the  two  SPS's  to  be  convolved  are  not  strictly  peri¬ 
odic,  but  are  actually  anti-periodic  depending  on  the  symmetric  extension  operator  for  that  case. 
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Table  2.  Inverse  Discrete  Trigonometric  Transforms 


Type 

Relation  to 

Forward  Matrix 

Input 

Output 

Equivalent 

Transform 

Symmetry 

Range 

Symmetry 

Range 

wsws 

O^N 

wsws 

^Ofi^WSiVS 

Su 

WAWA 

WAWA 

1  -^W-l 

WSWA 

O^N-l 

HSHS 

0-^W-l 

^qI^WSWA 

Sll 

WAWS 

\-^N 

HAHA 

0-^W-l 

'~J0Ol^WAWS 

HSHS 

O^N-l 

WSWA 

0-»W-l 

01fi^HSHS 

HAHA 

0->N-l 

WAWS 

1  ->N 

"^J^lfi^HAHA 

e;l 

HSHA 

0-»W-l 

HSHA 

0-^N-l 

^x\^HSHA 

SaI 

HAHS 

O^W-1 

HAHS 

O^N-l 

*“7^1  \^HAHS 

2’2 

e^o 

C-'-  *  c 

lN-\^ 

WSHS 

0-^N-l 

WSHS 

0-^N-l 

0o!o^wshs 

5,;' 

WAHA 

l^N-l 

WAHA 

^J^Ofi^WAHA 

^0 

c'=2;_i^3. 

WSHA 

0->N-l 

HSWS 

0->N-l 

0q1^WSHA 

si 

WAHS 

l^N-l 

HAWA 

O-^N-2 

1 

I 

Czl 

HSWS 

0->  A-1 

WSHA 

0-^N-l 

si 

HAWA 

0->N-2 

WAHS 

1  ->W-1 

^J^Xfi^HAWA 

C~^  -  ^  c 

IN -\*° 

HSWA 

O^N-2 

HSWA 

O^W-2 

^XX^HSWA 

si 

^  C 

HAWS 

0->W-l 

HAWS 

0^7V-1 

~j^x\^HAWS 

The  formula  for  circular  convolution,  denoted  by  ©,  is  well-known  [37].  Skew-circular  convo¬ 
lution  is  expressed  as 


d{n)  =  0{ri)(i)  h{ri) 

n 

=  d{m)h{n  -m)-  0(m)h(n  -m  +  N). 


N-l 


m=0 


m=n+\ 


Martucci's  key  result  [34]  is  that  Eq.  (3 1)  is  exactly  equivalent  to 
din  -  «o)  =  Z'  {7i  {Bin)}  •  7^  {/»(«)}}, 


(32) 


(33) 
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where  denotes  one  of  the  16  DTT  operators.  The  value  of  equals  either  0  or  1  indicat¬ 
ing  that  in  some  cases  the  multiplication  and  inverse  trigonometric  transform  operations  intro¬ 
duce  a  one-sample  delay  in  the  output.  Table  3  lists  all  forty  cases  from  [34]  of  the  symmetric 
convolution-multiplication  property,  and  identifies  those  cases  where  the  output  is  delayed  by 
one  sample. 


Table  3.  Forty  Cases  of  Symmetric  Convolution-Multiplication 


1  N  even,  M  =  IN 

1  Noddi,M  =  2N-l  \ 

Si 

Si 

IDl 

no 

Si 

Si 

@ 

"7, 

m 

IBB 

no 

Swsws 

Swsws 

© 

Cl. 

Si. 

0 

SwSHS 

SwSHS 

© 

Su 

e;l 

0 

Swsws 

SwAWA 

© 

Cl, 

Su 

Su 

0 

SwSHS 

SwAHA 

© 

Cl. 

5.. 

s:l 

0 

SwAWA 

SwAWA 

© 

Su 

Su 

-e:: 

0 

SwAHA 

SwAHA 

© 

^\o 

-e:l 

0 

Shshs 

Swsws 

© 

Su 

Su 

ei] 

0 

Shsws 

SwSHS 

© 

e^. 

Cu 

eil 

0 

Shshs 

SwAWA 

© 

e^. 

Su 

Su 

0 

Shsws 

SwAHA 

© 

e^. 

S\o 

Sll 

0 

Shaha 

Swsws 

© 

Si. 

Si, 

a 

0 

SwSHS 

© 

^2o 

ei„ 

BB 

0 

Shaha 

SwAWA 

© 

Su 

Su 

-e;] 

0 

Shaw  A 

SwAHA 

© 

Su 

-et 

0 

Shshs 

Shshs 

© 

e^, 

e^, 

1 

Shsws 

Shsws 

© 

eu 

e;l 

1 

Shshs 

© 

e^, 

S^, 

1 

Shsws 

Shaw  A 

© 

Su 

s;l 

1 

ISRI 

Shaha 

© 

Su 

^le 

1 

Shaw  A 

Shaw  A 

© 

Su 

Su 

-e;J 

1 

^ff'SfVA 

SwSWA 

© 

Su 

0 

SwSHA 

SwSHA 

© 

Cu 

Cu 

e;l 

0 

^WSfVA 

SwAWS 

© 

Su 

Bl 

0 

SwSHA 

SwAHS 

© 

Cu 

Su 

0 

SwAWS 

SwAWS 

© 

Su 

-e;l 

0 

SwAWS 

SwAHS 

© 

Su 

S-io 

0 

^HSHA 

SwSWA 

© 

Su 

e^. 

e;] 

0 

Shswa 

SwSHA 

© 

Cu 

Cu 

Su 

0 

Shsha 

mm 

© 

Cu 

Su 

0 

Shswa 

SwAHS 

© 

Su 

0 

Shahs 

SwSWA 

© 

Su 

Cu 

0 

Shaws 

SwSHA 

© 

^4o 

Cu 

0 

Shahs 

SwAWS 

© 

Su 

Su 

-e:l 

0 

Shaws 

SwAHS 

© 

•^4o 

Su 

-e;l 

0 

Shsha 

Shsha 

© 

Su 

Su 

1 

Shswa 

Shswa 

© 

Cu 

Cu 

e;l 

1 

Shsha 

Shahs 

© 

Su 

Su 

1 

Shswa 

Shaws 

© 

Su 

•^4o 

1 

Shahs 

Shahs 

© 

Su 

Su 

-e;l 

1 

Shaws 

Shaws 

© 

Si,  a 

54. 

-et 

1 

Each  group  of  transforms  in  an  entry  of  Table  3  performs  convolution  in  the  sequence  do¬ 
main  by  transforming  the  two  sequences,  multiplying  the  results  in  the  transform  domain,  and 
then  taking  the  inverse  transform  to  produce  an  answer  in  the  sequence  domain.  Using  the  sym- 
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metric  convolution-multiplication  property  results  in  a  savings  in  computational  complexity  as 
fast  algorithms  exist  for  the  DTTs  [5],  [39],  [50],  The  algorithms  have  the  same  computational 
complexity  as  the  fast  Fourier  transform. 

The  symmetric  convolution-multiplication  property  of  the  DTTs  requires  implicit  symme¬ 
try  in  the  sequences  being  convolved.  The  underlying  symmetry  is  analogous  to  the  implied  pe¬ 
riodicity  required  of  sequences  that  are  circularly  convolved  [37].  Martucci  states  [34]  that  the 
symmetric  convolution-multiplication  property  he  derives  for  one-dimensional  sequences  extends 
to  asymmetric  sequences  if  they  are  first  decomposed  into  their  symmetric  and  antisymmetric 
parts.  This  decomposition  is  straightforward  in  one  dimension,  but  can  be  complicated  in  multi¬ 
ple  dimensions  because  an  asymmetric  sequence  must  be  decomposed  into  its  symmetric  and  an¬ 
tisymmetric  parts  in  each  dimension  across  multiple  lines,  planes,  or  hyperplanes  of  symme¬ 
try  [34]. 

The  first  section  of  this  background  chapter  described  how  DFT  matrices  diagonaliae  cir¬ 
cular  convolution  matrices.  The  next  chapter  shows  how  to  derive  similar  relationships  where 
the  DTT  matrices  presented  in  this  last  section  diagonalize  symmetric  convolution  matrices. 
When  represented  by  matrices,  the  symmetric  convolution-multiplication  property  extends  more 
easily  to  both  multidimensional  and  asymmetric  sequences.  The  new  matrix  form  of  the  sym¬ 
metric  convolution-multiplication  property  then  becomes  the  basis  for  deriving  inverse  and  scalar 
Wiener  filters  in  the  trigonometric  transform  domain.  These  filters  are  shown  to  be  similar  but 
improved  versions  of  the  Fourier  domain  representations  presented  in  the  second  section  of  this 
chapter. 
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III.  Symmetric  Convolution  of  Asymmetric 
Multidimensional  Sequences 


An  alternate  method  of  deriving  Martucci's  [34]  symmetric  convolution-multiplication 
property  of  discrete  trigonometric  transforms  is  presented  in  this  chapter.  A  vector-matrix  form 
of  the  property  is  produced  by  this  alternate  method  [9].  A  new  procedure  is  demonstrated  of 
how  the  convolutional  forms  of  trigonometric  transforms  can  operate  as  matrices  which  diagonal¬ 
ize  symmetric  convolution  matrices  [11].  The  diagonalized  forms  of  the  symmetric  convolution 
matrices  then  allow  the  symmetric  convolution-multiplication  property  to  be  easily  extended  to 
multidimensional  and  asymmetric  sequences. 

3.1  Vector-Matrix  Derivation  of  the  Symmetric  Convolution-Multiplication  Property 

It  is  possible  to  derive  the  symmetric  convolution-multiplication  property  for  discrete 
trigonometric  transforms  (DTTs)  using  vector-matrix  methods  [9].  The  result  of  the  derivation  is 
very  similar  to  Hunt's  result  [23]  for  the  discrete  Fourier  transform  (DFT)  presented  as  back¬ 
ground  in  the  previous  chapter.  This  alternate  derivation  provides  a  new  look  at  the  symmetric 
convolution-multiplication  property  for  the  DTTs,  just  as  Hunt's  result  provided  a  vector-matrix 
proof  of  the  circular  convolution-multiplication  property  for  the  DFT.  The  derivation  which 
follows  relies  heavily  on  his  key  result  expressed  as  Eq.  (1). 

The  first  step  in  the  vector-matrix  derivation  of  the  symmetric  convolution-multiplication 
property  is  to  recognize  the  matrix  relationship  between  the  DFT  and  the  generalized  DFT 
(GDFT).  An  M  X  M  DFT  matrix,  W^,  defined  in  Eq.  (2)  is  expressible  in  terms  of  each  of  the 
four  forward  and  inverse  GDFT  matrices,  ^  where  a  and  b  equal 

0  or  y.  The  DFT  matrix  relates  to  each  of  these  GDFT  matrices  through  an  Mx  M  diagonal 
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matrix  with  elements  =exp{y^/M}  for  m  =  0,\,...,  M-l.  The  sizes  of  the 

transform  matrices  and  the  matrix,  Df^ ,  which  relates  them  are  all  M  x  M  since  they  operate  on 
symmetric  extensions  of  the  sequences  to  be  convolved.  Recall  that  the  original  sequences  to  be 
convolved  had  length  N,  but  the  symmetric  extension  required  to  perform  symmetric  convolution 
increases  their  length  to  M  =  2N  for  iV  even  or  M^2N-\  for  7/ odd.  The  relationship  among 
the  four  forward  GDFTs  is 


~  ^mGq  j 


The  inverse  DFT  matrix  is  expressible  in  terms  of  the  inverse  GDFT  matrix  for  each  case  by  in¬ 
verting  the  above  expressions.  All  inverses  exist  so  that 


~  Gq  o  m 


7)-*  -  p~j' 


An  additional  convolution-multiplication  property  exists  for  the  skew-circular  convolution 
operation  d  =  h®d  =  Hsd,  where  is  an  iV x  A  skew-circulant matrix.  That  is, 


h{<S)  -h{N-l)  •••  -/?(!) 

hi\)  hid)  •••  -h{2) 

h{N-\)  h(N-2)  ■■■  hiO) 


(36) 


As  mentioned  previously,  skew-circular  convolution  is  the  underlying  form  of  convolution  in 
half  of  the  forty  cases  of  symmetric  convolution  [34].  The  convolution-multiplication  property 
for  skew-circular  convolution  is  an  extension  of  a  result  of  Vemet's  [49].  His  result  begins  with 
the  GDFT  of  0(n)  for  ^  and  b  =  0, 


=  ^  0(n)  exp(^-72;r^^i^j , 


(37) 
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for  A:  =  0,  1,  . . - 1.  The  subscript  'O'  in  the  term  indicates  the  transform  domain  of 

the  odd  DFT  as  Vemet  defines  this  transform.  The  skew-circulant  convolution  of  two  sequences, 
h{n)  and  d{n),  is  exactly  equivalent  to  the  inverse  odd  DFT  of  the  product  of  their  transforms, 
or  equivalently. 


J_ 

N 


rv  -1  .  . 

^  ‘^o^k)9o{k)  exp(^y2;r^^i^^ 

A:=0 


(38) 


for  «  =  0,  1,  . . .,  A'"  - 1.  In  terms  of  matrices,  Martucci  [34]  generalizes  Vemet's  result  for  the 
inverse  odd  DFT  to  =i^o.j.Ar- 

Equation  (38)  implies  that  an  odd  DFT  matrix  diagonalizes  a  skew-circulant  matrix 
through  the  relation  •  The  odd  DFT  matrix ,  ,  equals  the  Nx  N  GDFT  ma¬ 
trix,  where  <3  =  ^  and  b  =  0.  The  matrix  is  diagonal  with  the  elements  of  G^Q^jh 

along  the  diagonal.  Thus,  for  skew-circular  convolution,  Gx^^^d  =  ^  ^ 


The  matrix  is  the  odd  DFT  matrix  which  operates  on  symmetrically  extended  se¬ 
quences  of  length  M  This  matrix  will  diagonalize  an  MxM  skew-circulant  convolution  ma¬ 
trix.  Like  the  DFT  matrix  in  Eqs.  (34)  and  (35),  the  matrix  is  expressible  in  terms  of  each  of 
the  four  special  cases  of  the  MxM  GDFT  for  which  a  and  b  equal  0  or  Here  the  relation¬ 
ships  are 


Again  all  inverses  exist,  so  that 


(39) 


(40) 
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The  matrices  in  Eqs.  (34)  and  (39)  and  their  inverses  in  Eqs.  (35)  and  (40)  will  pre-  and 
post-multiply  the  circulant  and  skew-circulant  matrices  underlying  symmetric  convolution  as  part 
of  deriving  the  symmetric  convolution-multiplication  property  of  the  DTTs.  The  result  of  this 
multiplication  must  be  diagonal  for  each  of  the  forty  cases  of  symmetric  convolution  to  show  that 
the  multiplication  property  holds.  First,  however,  symmetric  convolution  must  be  expressed  as  a 
matrix  operation.  Equation  (3 1)  which  used  operator  notation  is  equivalent  to 

d^R,H^E,e,  (41) 

which  is  expressed  using  vector-matrix  notation.  The  vectors  Oandd  are  Zxl  input  and  output 
vectors.  The  matrix  Ei  is  an  Mx  I  symmetric  extension  matrix.  The  matrix  is  an  MxM 
circulant  or  skew-circulant  matrix  with  the  vector  E^h  as  its  first  column.  The  matrix  = 

\li  O],  where  is  an  Lx  L  identity  matrix,  and  0  is  an  appropriately-sized  zero  matrix  to 
make  Rj^mLxM  matrix.  The  length  L  equals  N -\,N,  ox  N  +  \  depending  on  the  type  of 
symmetry  implied  and  the  different  sizes  for  input  and  output  sequences  listed  in  Table  1.  As 
before,  the  length  Mequals  2N  for  N  even  ox  2N -\  for  A  odd. 

All  of  the  symmetric  extensions,  shifts,  additions,  and  multiplications  of  symmetric  convo¬ 
lution  in  Eq.  (41)  can  be  captured  in  a  single  Lx  L  matrix,  =  Rj^HyE^ ,  so  that  d  =  H^9. 
The  subscript  'SC  denotes  symmetric  convolution,  which  is  generalized  in  this  derivation  to  any 
one  of  the  forty  cases  of  symmetric  convolution.  The  form  of  the  matrix  Hsc  will  be  different 
for  each  of  the  forty  different  cases.  The  goal  of  this  alternate  method  of  deriving  the  symmetric 
convolution-multiplication  property  is  to  show  that  an  LxL  diagonal  matrix  form,  exists 
for  the  symmetric  convolution  matrix,  Hsc-  The  subscript  T  denotes  that  the  diagonal  matrix 
exists  in  the  trigonometric  transform  domain.  It  will  be  shown  that  LxL  trigonometric 
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transform  matrices,  represented,  in  general,  by  the  matrix  cause  this  diagonalization.  The 
matrix  can  represent  any  one  of  the  sixteen  DTTs. 

Before  proceeding  with  the  derivation,  it  is  useful  to  present  a  summary  of  the  three  types 
of  one-dimensional  transforms  discussed  thus  far  in  tabular  form.  This  summary  for  the  DFT, 
the  odd  DFT,  and  an  arbitrary  DTT  appears  in  Table  4.  All  of  the  matrices  related  to  the  DFT 


Table  4.  One-Dimensional  Transforms 


BH 

Matrix  Repre¬ 
sentation 

Uncierlying  Form 
of  Convolution 

Convolution 

Matrix 

Diagonal  Form  in 
Transform  Domain 

DFT 

circular 

He 

odd-DFT 

I/"l 
^  M 

skew-circular 

Hs 

DTT 

symmetric 

Hsc 

and  the  odd  DFT  listed  in  Table  4  have  dimension  Mx  M  because  they  operate  within  the 
definition  of  symmetric  convolution  on  symmetric  extensions  of  the  sequences  being  convolved. 
The  trigonometric  transform  matrix,  T^,  and  the  matrices  Hgc  and  all  have  dimension 
LxL  because  the  individual  DTT  matrices  can  have  dimension  N  -\x  N  -I,  Nx  N,  or 
N  +  lxN  +  l  depending  on  the  choice  of  transform  from  Table  1  in  the  previous  chapter.  The 
dimension  of  each  of  the  sixteen  transforms  is  based  on  the  underlying  type  of  symmetry  present 
in  the  sequence  being  transformed.  Table  1  also  listed  the  underlying  symmetry  and  sequence 
length  for  each  of  the  sixteen  DTTs. 

« 

To  proceed  with  the  vector-matrix  derivation  of  the  symmetric  convolution-multiplication 
property,  the  diagonal  form  of  the  matrix  must  next  be  substituted  into  Eq.  (41).  The  matrix 
will  be  either  circulant  or  skew-circulant  depending  on  the  underlying  form  of  convolution 
in  the  particular  case  of  symmetric  convolution  from  Table  3  that  Eq.  (41)  represents.  Recall  that 
Eq.  (41)  generically  represents  any  one  of  the  forty  cases.  Thus  the  appropriate  substitution  for 
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in  Eq.  (41)  will  be  where  “^p  =  diag|fF^'JE'2A|  if  the  underlying  form  of  con¬ 
volution  is  circular.  The  substitution  for  H)^  in  Eq.  (41)  will  be  where  = 

Ai^g\v-jE,h]  if  the  underlying  form  of  convolution  is  skew-circular.  The  diag{»}  operation 


creates  a  diagonal  matrix  from  a  vector.  Any  of  the  equivalent  expressions  from  Eqs.  (34),  (35), 
(39),  or  (40)  then  become  allowable  substitutions  for  W^,  W^,  V^,  or  Vm  in  the  result.  With 
these  substitutions,  Eq.  (41)  becomes 

=  RpG;ymg{G^E,h]G,E,d]  (42) 

The  matrix  in  Eq.  (42)  is  either  “^p  or  depending  on  the  underlying  form  of  convolution. 
The  three  matrices,  each  represent  one  of  the  four  special  cases  of  an  Mx  M  GDFT  matrix 
for  which  a  and  b  equal  0  or  y.  In  some  cases  the  matrix  G„  may  include  multiplication  by  a 
factor  of  j  or  -  j.  This  factor  is  necessary  because  of  the  relationship  between  each  trigonomet¬ 
ric  transform  and  the  GDFT. 

The  relationship  between  each  DTT  operator,  7{»},  and  the  GDFT  operator, 

shown  in  Tables  1  and  2  for  the  forward  and  inverse  transforms.  One  of  Martucci’s  results  [34]  is 
that  the  DTT  of  a  sequence  can  be  expressed  as  a  GDFT  operating  on  a  symmetrically  extended 

version  of  the  sequence.  That  is,  7{x{n)]=  where  £se{»}  is  one  of  the  sixteen 

symmetric  extension  operators,  hence  bearing  the  subscript  'SE'.  The  example  was  given  in  the 
previous  chapter  that  <^,{x(«)}  =  ^o.o{%s3rs{^(”)}}-  Another  example  which  shows  the 
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equivalence  of  a  DTT  to  the  GDFT  is  S'4o{j:(«)}  from  Table  2.  Here  the 

need  for  multiplication  by  a  factor  of  -  j  is  illustrated. 

An  equivalence  between  DTT  and  GDFT  matrices  exists  that  is  similar  to  the  equivalence 
between  DTT  and  GDFT  operators.  If  the  sequence  x{n)  is  expressed  as  an  Z  x  1  vector  x,  then 

it  can  be  symmetrically  extended  by  the  Mx  L  matrix  ^SE  as  demonstrated  in  the  previous 
chapter.  This  extension  produces  the  Mxl  vector  x  =  E^eX.  Now  multiplying  the  vector  x 
by  an  Mx  M  GDFT  matrix,  *  a/ »  will  be  exactly  equivalent  to  multiplying  the  vector  x  by 

the  MxL  extended  DTT  matrix  f.  That  is,  fx  =  Ggj,  ,^x  =  G„  ,,  i^EeeX,  so  that  the  vector- 
matrix  equivalent  form  of  7{x(«)}  =  ^aj){^SE  {^(”)}}  tilde,  ~,  on  the 

transform  matrix,  f,  indicates  that  it  is  an  extended-length  transform.  Normally  a  DTT  matrix, 
Te,  has  dimension  LxL,  but  its  Mx  L  extended  version,  T,  has  twice  as  many  rows.  The 
equivalence  between  DTT  and  GDFT  matrices  allows  substitutions  of  the  form  T„  =  G„E^,  to  be 
made  in  Eq.  (42).  These  substitutions  produce 

d  =  REG^^^fjh  ®  (43) 

The  portion  of  Eq.  (43)  inside  the  brackets  results  in  an  M  x  1  vector  because  of  the  point-wise 
multiplication.  The  vectors  dandh  both  have  dimension  Zxl,  but  are  transformed  into  Mxl 

vectors  by  the  MxL  transformation  matrices  and  T^.  The  Mx  1  vectors  T^d  and  will 
both  represent  symmetric  periodic  sequences  (SPS's).  They  will  be  SPS's  because  the  DTTs  re¬ 
quire  an  implied  symmetry  associated  with  both  their  input  and  output  sequences.  Just  like  the 
input  and  output  sequences  to  a  DFT  are  periodic,  the  input  and  output  sequences  to  a  DTT  are 
symmetric  periodic.  The  product  of  any  two  SPS's  is  itself  an  SPS,  so  that  the  point-wise  product 
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fjh  ©  f^O  will  also  represent  an  SPS.  Because  the  product  fjh  0  f^O  has  dimension  Mx\,  it 
will  have  twice  as  many  samples  as  the  L  samples  which  are  needed  to  completely  characterize  it 
as  an  SPS.  The  L  samples  which  completely  characterize  the  vector  ©  f^O  are  analogous  to 
the  samples  which  characterize  a  periodic  sequence  by  its  fundamental  period.  Since  the  vector 
T2h  ©  T^O  can  be  completely  characterized  by  Z  samples,  the  MxL  matrices  7J  and  2^  in 
Eq.  (43)  can  be  reduced  to  their  standard  Lx  L  size  represented  by  and  T2 .  A  third  Mx  L 
symmetric  extension  matrix,  E^,  can  then  symmetrically  extend  the  result  of  the  Z  x  1  point- 
wise  multiplication  T2h  ©  T^O.  Making  the  substitution  T^A  ©  =  E^\T2h  ©  in 

Eq.  (43)  produces 


d  =  RLG;^E2[T2h  ©  T^e] 
=  RJ{\T2h  ©  T,e], 


(44) 


by  recognizing  that  G^^E^  =  Zj”*.  The  effect  of  the  Lx  M  matrix  is  to  retain  the  first  Z 
samples  of  the  Mx  1  vector  2^“*[Z2A  ©  Z,^].  These  Z  samples  are  the  same  which  result  from 
reducing  the  Mx  Z  matrix  in  Eq.  (44)  to  its  standard  Lx  L  size  represented  by  Z3.  Thus 


rf  =  Z3"'[Z2A  ©Z,4  (45) 

where  2],  7^,  and  Zj"'  represent  Lx  L  trigonometric  forward  and  inverse  transforms.  Equa¬ 
tion  (45)  is  the  symmetric  convolution-multiplication  property  of  the  discrete  trigonometric 
transforms  expressed  in  vector-matrix  notation.  An  alternate  and  sometimes  more  useful  method 
of  expressing  the  result  in  Eq.  (45)  is 

d  =  T;^‘»jT^e.  (46) 
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The  Z,  X  Z  diagonal  matrix  equals  diagjT^Zr},  where  the  subscript  T'  again  indicates  the 
trigonometric  transform  domain. 

The  preceding  general  procedure  allows  for  the  derivation  of  any  one  of  the  forty  versions 
of  symmetric  convolution  in  Table  3  using  vector  notation.  The  key  to  deriving  any  particular 
case  is  to  make  the  appropriate  substitutions  for  each  and  <7„  in  Eqs.  (41)  -  (44)  which  will 

yield  the  correct  three  transforms,  7^,  in  Eqs.  (45)  and  (46).  The  appropriate  substitutions  for 
and  for  each  of  the  forty  cases  appear  in  Table  5.  Note  that  the  entries  in  the  table  tem¬ 
porarily  suppress  the  subscript  'M'  indicating  the  dimension  of  the  transform  since  all  GDFT  en¬ 
tries  have  dimension  My.  M.  Note  that  some  inverse  transforms  in  the  table  have  the  form 

^0 1 M  M’  where  the  T'  indicates  a  whole-sample  delay.  These  delaying  forms  of  the 

GDFT  result  from  the  expressions  m  ~  i  m 

each  case  where  a  one-sample  delay  occurs  in  the  symmetric  convolution  output  in  Eq.  (33)  for 
the  cases  where  =  1. 

As  an  example  of  how  to  derive  one  of  the  forty  types  of  symmetric  convolution,  consider 
the  case  where  the  N  y  \  vectors  h  and  6  have  WSWS  and  HSHS  symmetry,  respectively.  In 
this  case 

d  =  R,HcE„SHsO,  (47) 

where  is  an  Mx  M  circulant  matrix  with  the  M x  1  vector  E^r^ws^  its  first  column.  The 

matrix  He  is  then  replaced  by  where  =  A\dig\w;}E^swsf^-  The  facts  from 

Eqs.  (34)  and  (35)  that  and  allow  He  to  be 
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Ewsws 

Ewsws 


Ewawa 

Ehshs 

Ehshs 

Ehaha 

\Ehaha 

Ehshs 

Ehshs 


Ehaha 

Ewswa 


Ewswa 

Ewaws 


Ehsha 

Ehsha 

\Ehahs 


Ehahs 


Ehsha 


Ehsha 


\Ehahs 


rewritten  as  Hq  —  where  Equation  (47)  then 


Table  5.  Substitutions  in  Eqs.  (41)  -  (44)  used  to 
Derive  the  Forty  Cases  of  Symmetric  Convolution 


N  even,  M=2N 


Nodd,M=2N-\ 


G 


g;' 


■'WSWS 


^wsws 


0,0 


WSHS 


■'WSHS 


■'WSHS 


G~^ 

0.0 


'WAWA 


-WAWA 


G. 


•/^O.O 


-jGZ 


'WSHS 


-WAHA 


'WAHA 


-WAWA 


'WSWS 


J^O.O 


j^o.o 


-G“' 

0,0 


'WAHA 


■'WAHA 


-WSHS 


0,0 

j^o.o 


JG. 


0,0 


-JK, 


-g: 


■'WSWS 


'WSWA 


G 


o.i 


-HSWS 


-WSHS 


-WSHA 


G 


0,4- 


0,0 


E, 


WAWA 


'WAWS 


o.i 


-HSWS 


-WAHA 


'WAHS 


o,i 


JG, 


0,0 


-jg:\ 


'WSWS 


-WAWS 


-jg: 


O.i 


-HAWA 


'WSHS 


-WAHS 


jG.i 


'WAWA 


'WSWA 


j^oA 


jGo,o 


-^oA 


'HAWA 


-WAHA 


-WSHA 


JG. 


oA 


JGo,o 


-^oA 


'WSWS 


0,4^ 


0.4 


g: 


'HSWS 


-HSWS 


'WSHS 


G 


0.4 


0,4 


G'^ 

0,1 


-HAHA 


'WAWA 


0,4 


-JGZ 


'HSWS 


-HAWA 


'WAHA 


0,4 


•^’^0.4 


-HAHA 


'WSWS 


JG.l 


JG, 


OA 


E. 


HAWA 


\Eh 


'WSHS 


JGoi, 


JG. 


o.i 


-JGZ 


‘WSWA 


‘HSHS 


G, 


i.o 


‘WSHA 


‘WSHA 


‘HSWS 


G, 


i.o 


G,. 


'WAWS 


'HAHA 


^40 


JG±o 


-K 


'WSHA 


'WAHS 


-HAWA 


% 


JGi. 


-K 


'WAWS 


'HSHS 


Me 


JG,, 


'WAWS 


'WAHS 


'HSWS 


J^Ao 


JGxo 


'WSWA 


'HSHA 


'WSHA 


'HSWA 


^,0 


'WAWS 


'HAHS 


G. 


J^Ao 


-jg; 


'HSWA 


'WAHS 


'HAWS 


JG,. 


■jg; 


'WSWA 


'HAHS 


Ml 


%0 


-K 


'HAWS 


'WSHA 


'HAWS 


JG„ 


i.o 


-JG'A 


'WAWS 


'HSHA 


Ml 


J^i.o 


-g:\ 


'HAWS 


'WAHS 


'HSWA 


JG„ 


JG,, 


-g;[ 


'HSHA 


'HSHS 


'HSWA 


'HSWA 


'HSWS 


^.1 


'HAHS 


'HAHS 


JG,, 


-J^t 


'HSWA 


'HAWS 


E^ 


HAWA 


Mi 


-K 


'HAHS 


'HSHS 


Ml 


JG,, 


-K 


'HAWS 


'HAWS 


'HSWS 


Ml 


JG„ 


becomes 


diag(  (^ofiM^wsws^  ^(>,\,m^hshs 

= 

[CuM^  ®  C2eM0] 

II 

^WSW^^le,N^  ® 

^\e,N^  ®  ^2e,N^\ 

“  ^2e,N^\e,N^  ®  ^le,N^\  (48) 

The  matrix  Rj^  retains  the  first // samples  of  the  Mx  1  vector  C2”Jm[Ci,  ©  Cj,  f]6^.  An  al¬ 
ternate  way  to  express  the  result  of  Eq.  (48)  is  in  the  matrix  form 

^  ~  ^2e,N^T,s^2e,N^>  (49) 

where  the  x  A  diagonal  matrix  =  diag|C]^  The  lowercase  subscript  's'  indicates  the 
vector  has  whole-sample  symmetry  in  the  transform  domain.  The  end  results  of  Eqs.  (48) 

or  (49)  are  exactly  equivalent  to  the  same  case  of  symmetric  convolution  expressed  using  opera¬ 
tor  notation  in  Table  3. 

By  making  the  appropriate  substitutions  from  Table  5,  the  remaining  39  cases  of  symmet¬ 
ric  convolution  are  just  as  easily  derived.  The  vector-matrix  derivation  presented  here  is  there¬ 
fore  equivalent  to  all  forty  cases  of  the  symmetric  convolution-multiplication  property  of  the 
DTTs  which  Martucci  [34]  derives  using  operator  notation. 

3.2  Extension  to  Multiple  Dimensions 

An  advantage  to  deriving  the  symmetric  convolution-multiplication  property  in  terms  of 
vectors  and  matrices  is  that  for  orthogonally-sampled  data  in  multiple  dimensions,  the  results  of 
the  previous  section  extend  naturally  using  Kronecker  products  as  defined  in  Eq.  (5).  Data  which 
are  orthogonally  sampled  in  two  dimensions  have  their  samples  aligned  on  a  unit-square  grid. 
Before  extending  the  results  of  the  previous  section,  diagonalizing  forms  for  the  two-dimensional 
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odd  DFT  are  presented.  The  two-dimensional  odd  DFT  is  used  in  cases  where  skew-circular 
convolution  is  the  underlying  form  of  convolution  in  symmetric  convolution.  These  diagonaliz¬ 
ing  forms  lay  the  groundwork  for  the  presentation  which  follows  on  two-dimensional  forms  of 
the  previous  section's  results.  The  case  of  multiple  dimensions  higher  than  two  is  also  consid¬ 
ered  in  this  section. 

A  result  similar  to  the  diagonalizing  form  of  two-dimensional  circular  convolution  in 
Eq.  (7)  exists  for  two-dimensional  skew-circular  convolution  based  on  the  two-dimensional  odd 
DFT.  In  this  case  Hg^  is  an  N^N2  block  skew-circulant  matrix  of  partitioned  blocks  oi 

skew-circulant  matrices  arranged  in  a  skew-circulant  pattern.  Odd  DFT  matrices  will  diagonalize 
the  two-dimensional  skew-circulant  convolution  matrix,  Hgg,  through  the  relation 

Hbs  =  {Vn,  ®  ®  VnI)-  (50) 

The  matrix  equals  diagj^F^^' where  h  =  vec^H^,  the  lexicographic  representation 

of  a  system  point  spread  function  (PSF),  H.  Thus  the  two-dimensional  skew-circular  convolu¬ 
tion-multiplication  property  is 

©(f^.'of;;)^,  (51) 

which  is  similar  to  Eq.  (8). 

The  possibilities  also  exist  for  a  matrix  of  partitioned  blocks  of  circulant  matrices  arranged 
in  a  skew-circulant  pattern,  or  a  matrix  of  partitioned  blocks  of  skew-circulant  matrices 

arranged  in  a  circulant  pattern,  The  #,1^2  x  matrices  Hg^c  Hc^g  operate  on 

lexicographically-ordered  vectors  which  represent  matrices.  The  subscript  'SxC  '  indicates  that 
the  matrix  operates  on  the  samples  corresponding  to  the  rows  using  circular  convolution 
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and  on  the  samples  corresponding  to  the  columns  using  skew-circular  convolution.  The  converse 
is  true  for  the  matrix  The  diagonalizing  forms  of  these  particular  cases  respectively  are 

Hs.c  =  K  <2)  ®  (52) 

and  (53) 

In  Eqs.  (52)  and  (53),  the  N^N2  x?/,A^2  matrices  and  ^p^o  ®qnal  diagj^F);^’ (8) 

and  diagj^IF;;^^*  respectively. 

The  problem  of  calculating  the  two-dimensional  symmetric  convolution  of  two  LiX  L2 
matrices,  H  and  0,  can  be  made  less  computationally  intense  by  the  above  diagonalizing  forms 
of  block  circulant  matrices,  block  skew-circulant  matrices,  or  combinations  of  the  two.  The  ma¬ 
trices  H  and  0  have  dimension  L^x  because  in  general  Z,,  equals  either  A^i  - 1 ,  A^i ,  or 
A^i  + 1  and  equals  either  N2  - 1,  N2,  or  N2  +  \  based  on  the  type  of  symmetric  extension 
underlying  the  rows  and  columns  of  H  and  0. 

An  equivalent  way  to  represent  the  matrices  Zf  and  0,  is  by  the  I^I^xX  lexicographi¬ 
cally-ordered  vectors  h  and  9.  The  two-dimensional  form  of  Eq.  (41)  is  thus 

(54) 

The  matrices  Ri^^  and  Rp^  have  the  same  form  as  from  before,  which  is  R^^^  =  oj  and 

Ri^  =[7^^  oj.  The  Z1Z2  ^  ^1^2  windowing  matrix  thus  retains  Z,  samples  from 

each  column  and  Z2  samples  from  each  row.  The  symmetric  extension  matrices  and  E^^ 
operate  on  the  rows  and  columns  of  0  separately  giving  it  possibly  different  symmetry  in  hori¬ 
zontal  and  vertical  directions.  The  dimensions  of  E^^  and  E^^  will  be  M,  x  Zj  and  M2  x  Z2,  re- 
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spectively,  where  M,  =  2N^  for  A^i  even  or  M,  =  2N^  - 1  for  AT,  odd,  and  M2  =  2N2  for  N2 
even  or  M2  =2A^2“1  for  ^2  odd.  The  xMiMj  matrix  Hg  can  be  either  block  circulant, 
block  skew-circulant,  or  a  combination  depending  on  the  underlying  symmetry  of  the  rows  and 
columns  of  the  L^xl^  matrix  H.  The  matrix  Hg  will  have  the  vector  {£2^  ®  its  first 

column.  The  matrices  Eic  and  will  also  have  dimension  Af,  x  Zj  and  M2X  L2,  respectively. 
They  will  perform  a  two-dimensional  symmetric  extension  on  the  matrix  H,  just  as  the  matrices 
and  perform  a  two-dimensional  symmetric  extension  on  the  matrix  0. 

As  before  in  the  one-dimensional  case,  all  of  the  symmetric  extensions,  shifts,  additions, 
and  multiplications  of  two-dimensional  symmetric  convolution  in  Eq.  (54)  can  be  contained  in  a 

single  L1L2  X  matrix  Hggc  =  (ifz,  ®  -^4  )i^b{E\c  ®  so  that  d  =  HggcO.  The  subscript 

'55'C' denotes  block  symmetric  convolution.  It  also  follows  from  the  one-dimensional  case  that  a 
diagonal  form,  exists  for  the  matrix  Hgg(..  The  diagonalization  in  two  dimensions  is  car¬ 
ried  out  using  Kronecker  products  of  the  DTTs  which  have  the  form  {t^  The  1^x1^ 

DTT  matrix  acts  on  the  samples  of  a  lexicographic  vector  which  represent  the  columns  of  a 
matrix.  Similarly,  the  Zj  x  Z2  DTT  matrix  T^  g^  acts  on  the  samples  which  represent  the  rows. 

A  summary  of  the  two-dimensional  transforms  discussed  thus  far  appears  in  Table  6.  All 
of  the  matrices  in  Table  6  related  to  the  DFT  or  odd  DFT  have  dimension  Mj  M2  x  M1M2  be¬ 
cause  they  act  on  symmetric  extensions  of  matrices  inside  the  definition  of  symmetric  convolu¬ 
tion.  The  two-dimensional  trigonometric  transform  matrix,  {T^g^®T^g^,  has  dimension 

Z,Z2  X  Z1Z2  because  it  acts  directly  on  a  lexicographic  vector  representing  a  matrix.  The  matri¬ 
ces  HggQ  and  in  Table  6  also  have  dimension  L^L^x  LJ-^. 
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The  procedure  to  extend  symmetric  convolution  to  multiple  dimensions  follows  a  similar 
derivation  as  the  previous  section's  derivation  for  one  dimension.  Here,  substitutions  from 
Eqs.  (34),  (35),  (39),  and  (40)  are  made  into  the  appropriate  diagonal  form  in  Eqs.  (6),  (50),  (52), 
or  (53).  The  DTT  matrices  which  result  from  these  substitutions  appear  in  the  two-dimensional 
equivalent  of  Eq.  (45)  as 

d  =  (55) 

An  alternate  way  to  express  Eq.  (55)  is 

D  =  0  (7;,07]^)]r3-/.  (56) 

The  two  results  in  Eqs.  (55)  and  (56)  are  equivalent  because  for  any  separable  transform  repre¬ 
sented  by  the  matrix  U,  the  vector  %  =  ®U^^x  is  lexicographically  equivalent  to  the  ma¬ 
trix  “Xu  =  [27].  The  DFT,  the  odd  DFT,  and  the  entire  family  of  trigonometric  trans¬ 

forms  are  all  separable.  In  Eqs.  (55)  and  (56)  the  set  of  DTTs,  7^,.,  which  acts  on  the 

rows  and  the  set  of  DTTs,  which  acts  on  the  columns  must  each  come  from  one 

of  the  allowable  sets  of  forty  transforms  in  Table  3,  but  they  need  not  be  the  same  set.  They  may 
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be  different  because  the  matrices  H  and  ©  may  have  different  underlying  symmetry  in  each  di¬ 
rection. 

These  same  diagonalizing  principles  apply  to  higher  dimensional  data  by  cascading  the 
Kronecker  products  in  the  diagonalizing  equations.  For  example  a  D-dimensional  circulant  con¬ 
volution  matrix  acting  on  D-dimensional  data  of  size  N^  x  is  diagonalized  by 

Hc.c.-.c=(Wn®W^®  -  -  OfF;;),  (57) 

where  is  diagonal  with  the  vector  ®  (E>  •••  along  the  diagonal.  In  this 

D-dimensional  case,  there  are  2^  possible  combinations  of  circulant  and  skew-circulant  multiple 
sub-blocks  in  the  block  circulant  structure  of  Up.  Also  in  this  multiple  dimensional  case,  the 
points  of  symmetry  underlying  the  sequences  in  one  dimension  and  the  lines  of  symmetry  under¬ 
lying  the  sequences  in  two  dimensions  become  planes  of  symmetry  in  three  dimensions  and  hy¬ 
perplanes  of  symmetry  in  more  than  three  dimensions.  In  the  next  section,  the  behavior  of  sym¬ 
metric  convolution  is  examined  for  cases  where  the  sequences  to  be  convolved  do  not  possess  the 
underlying  symmetry  required  for  the  symmetric  convolution-multiplication  property  of  the 
DTTs  to  hold. 

3.3  Extension  to  Asymmetric  Sequences 

Even  though  the  symmetric  convolution-multiplication  property  of  the  DTTs  extends  eas¬ 
ily  to  multiple  dimensions,  it  is  still  limited  by  the  underlying  symmetry  in  the  sequences  to  be 
convolved  or  multiplied  in  the  transform  domain.  For  multidimensional  sequences,  there  must  be 
symmetry  in  every  dimension  [34].  Martucci  mentions  [34]  that  asymmetric  sequences  can  be 
decomposed  into  their  symmetric  and  antisymmetric  parts,  transformed  using  different  trans¬ 
forms  because  of  the  different  underlying  symmetry,  and  then  multiplied. 
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This  decomposition  is  straightforward  for  one-dimension,  but  can  be  complicated  for  two 
or  more  dimensions.  A  finite  one-dimensional  sequence,  x{n),  which  is  zero  outside  the  region 
0  ^  «  <  - 1,  can  be  decomposed  as  x{n)  =  x^{n)  +  x^n)  in  two  ways.  Using  whole-sample 

symmetry. 


Xain)  = 


-\x{-n), 

0, 

\x{n). 


n<0 

^x(-n). 

n  <0 

n  =  0 

and  x^(n)  =  < 

x(n). 

n  =  0 

(58) 

«  >  0, 

jx(n). 

n>0. 

where  x^{n)  has  WA  symmetry  to  the  left  and  x^n)  has  WS  symmetry  to  the  left.  The  symme¬ 
try  to  the  right  for  x^in)  can  be  any  one  of  the  four  types  (WA,  WS,  HA,  or  HS)  as  long  as  the 
symmetric  extension  of  x^n)  to  the  right  cancels  with  it  to  yield  x{n)  =  x^(n)  +  x^(n)  which  is 
zero  outside  the  region  0  <  n  <  W  - 1.  For  example  if  x^(n)  has  WAHS  symmetry,  then  x^(n) 
must  have  WSHA  symmetry.  Similarly  using  half-sample  symmetry  to  the  left. 


,  .  \-2x{-n-\),  «<0  \\x{-n-\),  «<0 

n  >0,  «>0. 


(59) 


Again  the  right-hand  symmetry  is  arbitrary  as  long  as  the  two  types  of  right-hand  symmetiy  for 
x„{n)  and  x^{n)  cancel  with  each  other  for  the  two  sequences. 

In  two  dimensions,  there  are  four  ways  to  decompose  an  orthogonally-sampled  sequence 
using  combinations  of  half-sample  and  whole-sample  symmetry.  Each  decomposition  must  have 
four  terms  consisting  of  symmetric  and  antisymmetric  parts  in  each  dimension  so  that 


x(«i,«2)  =  X^„(«1,«2)  +  X^,(«i,n2)  + 


(60) 


Definitions  of  the  four  terms  for  each  of  the  four  possible  decompositions  appear  in  Table  7.  The 
orthogonal  coordinate  system  axes  form  the  lines  of  symmetry  because  of  orthogonal  sampling. 
These  definitions  exist  for  a  sequence  in  the  first  quadrant,  i.e.,  0  <  <  W,  - 1,  0  <  «2  ^  - 1. 
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Table  7.  Decomposition  of  a  Two-Dimensional 
Asymmetric  Sequence  into  Symmetric  and  Antisymmetric  Parts 


Similar  definitions  exist  for  sequences  in  other  quadrants.  This  decomposition  still  applies  to 
sequences  which  exist  in  more  than  one  quadrant.  The  decomposition  of  a  multiquadrant  se- 
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quence  requires  it  to  be  separated  into  subsequences  which  exist  in  the  individual  quadrants  of 
the  -  «2  plane.  The  four  subsequences  must  then  be  individually  decomposed  into  their  anti¬ 
symmetric-antisymmetric,  antisymmetric-symmetric,  symmetric-antisymmetric  and  symmetric- 
symmetric  parts,  and  the  results  added. 

These  concepts  apply  to  orthogonally-sampled  asymmetric  sequences  in  three  and  higher 
dimensions  as  well.  A  three-dimensional  sequence  must  be  decomposed  into  eight  symmetric 
and  antisymmetric  parts  in  each  dimension.  In  general  a  Z)-dimensional  sequence  will  have  2^ 
components  in  its  decomposition.  A  three-dimensional  sequence  will  have  planes  of  symmetry, 
and  a  Z)-dimensional  sequence  will  have  hyperplanes  of  symmetiy  [34].  This  multidimensional 
decomposition  is  very  similar  to  the  decomposition  needed  to  implement  a  multidimensional  Hil¬ 
bert  transform  [19],  except  here  the  phase  is  reversed  by  180°  instead  of  being  rotated  by  90°. 

In  the  previous  discussion  an  alternate  way  to  derive  Martucci's  symmetric  convolution- 
multiplication  property  [34]  was  shown  using  vector-matrix  notation  and  the  property  was  ex¬ 
tended  to  the  more  general  class  of  asymmetric  multidimensional  sequences.  The  focus  now 
turns  to  filtering  asymmetric  multidimensional  sequences  in  the  trigonometric  transform  domain. 

3.4  A  Filtering  Example 

A  problem  that  requires  the  convolution  of  multidimensional  asymmetric  sequences  is 
modeling  the  effects  of  atmospheric  turbulence  on  imaging  systems  [40].  This  model  provides  a 
good  example  of  a  problem  which  requires  multidimensional  asymmetry  to  demonstrate  the 
symmetric  convolution-multiplication  property  for  DTTs  because  models  for  atmospheric  turbu¬ 
lence  result  in  asymmetric  point  spread  functions  (PSFs). 

Figure  3  shows  an  example  of  a  256  x  256  PSF  selected  to  model  the  effects  of  atmos¬ 
pheric  turbulence  on  the  aperture  of  an  imaging  system.  The  entire  PSF  appears  in  Figure  3(a). 
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(e)  Decomposed  PSF,  central  samples, 
symmetric,  «2  antisymmetric 


(f)  Decomposed  PSF,  central  samples, 
symmetric,  «2  symmetric 


Figure  3.  Decomposition  of  Point  Spread  Function  (PSF) 
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The  central  32  x  32  samples  appear  in  Figure  3(b)  where  the  asymmetry  is  clearly  evident.  A 
direct  transformation  of  this  PSF  into  the  trigonometric  transform  domain  is  not  possible  because 
of  its  asjmimetry.  The  trigonometric  transforms  can  only  act  directly  on  the  decomposed  S5mi- 
metric  and  antisymmetric  parts  of  the  PSF  about  both  the  w,  and  «2  axes.  If  /?(«,,  n2)  denotes  the 
PSF,  the  decomposition  becomes 

/j(«l ,  «2  )  =  /Jaa  («i ,  «2  )  +  («1 ,  «2  )  +  («1  >  «2  )  +  («1 .  «2  )•  (61) 

The  decomposition  in  Eq.  (61)  is  not  direct  since  A(/7i,«2)  exists  in  all  four  quadrants  of  the 
-  «2  plane.  The  correct  decomposition  technique  is  to  decompose  the  portion  of  the  PSF  in 
each  quadrant  separately  and  then  add  the  results  having  like  symmetry.  For  example,  the  term 
in  Eq.  (61)  which  is  antisymmetric  in  «,  and  symmetric  in  «2  would  arise  from  the 
sum  of  the  antisymmetric-symmetric  portions  of  the  whole-sample  symmetric  decompositions  for 
each  quadrant.  The  central  samples  of  the  decomposition  of  this  PSF  using  whole-sample  sym¬ 
metry  for  both  «i  and  02  appear  in  Figures  3(c)-(f).  The  antisymmetric-antisymmetric  portion  of 

the  decomposition,  h^(rii,n2),  appears  in  Figure  3(c),  /i^(nj,«2)  appears  in  Figure  3(d), 
/)'^^(«,,«2)  appears  in  Figure  3(e),  and  h^^(n^,n2)  appears  in  Figure  3(f).  Half-sample  symmetry 

for  either  or  both  directions  could  have  just  as  easily  been  chosen  for  this  decomposition  instead 
of  whole-sample  symmetry.  The  only  impact  to  this  example  would  be  the  choice  of  trigono¬ 
metric  transform  to  apply  later. 

Figure  4  shows  the  negative  of  a  256  x  256  pixel  computer-generated  rendering  of  an 
ocean  reconnaissance  satellite.  The  goal  of  this  filtering  example  is  to  convolve  this  satellite 
object  with  the  PSF  in  Figure  3  by  converting  each  to  the  trigonometric  transform  domain,  point- 
multiplying  the  results  and  then  inverse  transforming  the  product  back  to  the  spatial  domain. 
Decomposing  the  object  to  be  imaged  into  its  symmetric  and  antisymmetric  parts  is  unnecessary 
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Figure  4.  Simulated  Satellite  Object 
(Negative  Shown  for  Clarity) 


because  the  convolution  shift  property  [37]  allows  for  a  shift  of  the  object's  origin.  The  convo¬ 
lution  shift  property  states  that  if  d^(n)  =  h{n)*0{n)  and  d2(n)=  h(n)* 0(n - Hq),  then  d2(n)  = 
-«q).  The  shift  of  the  object's  origin  allows  the  entire  object  to  appear  as  if  it  were  one 
fundamental  period  of  a  two-dimensional  symmetric  periodic  sequence.  The  symmetric  periodic 
extensions  of  the  object  are  external  to  the  object  itself  and  therefore  transparent  to  the  problem. 
The  symmetric  periodic  extensions  of  the  PSF  must,  however,  occur  internally  within  the  subse¬ 
quences  and  thus  need  to  be  accounted  for. 

The  DTTs  imply  an  underlying  symmetric  periodicity  in  both  the  sequence  and  the  trans¬ 
form  domains  in  the  same  way  that  the  DFT  of  a  finite  sequence  implies  periodicity  in  both  do¬ 
mains.  Therefore  the  implementation  of  a  filter  in  the  transform  domain  does  not  require  the  full 
symmetric  extensions  of  either  the  object  or  the  PSF.  The  symmetry  is  implied  by  the  trigono¬ 
metric  transforms  just  like  the  DFT  implies  periodicity.  In  fact  since  the  symmetry  of  the  four 
parts  of  the  decomposition  of  the  PSF  must  be  accounted  for  internally  for  each  subsequence,  the 
filter  needs  only  to  retain  the  principal  128  x  128  values  in  the  first  quadrant  of  each  part.  These 
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parts  in  the  first  quadrant  are  designated  by  A" («,, 772),  h"{n^,n2),  4o(«i,«2)>  and  AJ^(n,,«2), 
where  the  superscript  'rr'  indicates  the  filter  right-half  terms  about  both  and  Uj  •  The  128  x 

128  right-half  subsequences  have  the  remainder  of  each  array  padded  with  zeros  so  their  full 
sizes  are  all  256  x  256.  This  process  is  the  two-dimensional  equivalent  of  retaining  just  the 
right-half  samples  of  the  filter  impulse  response  and  then  zero-padding  in  the  one-dimensional 
examples  of  [34]  and  [35]. 

With  all  underlying  symmetry  properly  accounted  for,  trigonometric  transforms  must  next 
be  applied  to  the  object  and  the  decomposed  PSF  to  perform  symmetric  convolution  via  multipli¬ 
cation  in  the  transform  domain.  In  this  example  the  finite  length  of  the  object  being  imaged  is 
even  in  both  directions,  i.e.,  Af,  =  A^2  =  ^56,  so  any  one  of  the  eight  even-length  DTTs  can  be 
applied  for  each  direction,  n-y  and  nj  ■  The  same  transform  does  not  need  to  be  applied  for  each 
direction.  For  this  example  let  the  matrix  ©  represent  the  satellite  object,  and  then  apply  a  type- 
n  DCT  to  both  the  rows  and  columns  of  (9.  The  transform  domain  representation  of  the  matrix 
©  thus  becomes  matrices  in  this  expression  have  dimension  256  x  256. 

Thus  far  the  filtering  problem  has  allowed  freedom  to  choose  any  type  of  symmetry  in  the 
decomposition  of  h{n^,n2).  There  has  also  been  freedom  to  choose  any  of  the  even-length  trans¬ 
forms  to  apply  to  the  object  being  filtered.  In  this  example,  whole-sample  symmetry  has  been 
chosen  in  each  direction  for  the  decomposition  of  A(ni,«2),  and  a  type-II  DCT  has  been  chosen 
to  transform  both  the  rows  and  columns  of  the  object.  These  two  choices  now  dictate  the  type  of 
transforms  to  apply  to  the  individual  components  of  the  decomposition  of  the  PSF,  h{n^,n2). 

The  choices  also  restrict  which  inverse  transforms  can  be  used  to  produce  the  convolved  image. 
From  Table  3  which  specifies  the  types  of  symmetric  convolution,  there  are  only  two  allowable 
even-length  cases  of  symmetric  convolution  which  have  a  type-II  DCT  operating  on  one  se- 
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quence  and  whole-sample  symmetry  in  the  other  sequence.  In  vector-matrix  form,  the  allowable 
cases  are  for  A  whole-sample  symmetric,  and  ®^2e,N^\ 

for  h  whole-sample  antisymmetric.  Thus  the  filter  needs  to  apply  a  type-I  DCT  and  a  type-I  DST 
to  the  appropriate  symmetric  and  antisymmetric  rows  and  columns  of  the  decomposed  subse¬ 
quences  of  h{n„n,).  If  HZ,  HZ,  HZ,  and  HZ  represent  the  128  x  128  principal  values  of 

zero-paddcd  to  dimcnsion  256x256,  then  the 

transforms  to  apply  are 


le.A'i" 


rrrr^T 


^T,as  “ 


,HZsL 


^T,sa  ~  ^le,Ni^Z^le,N2  ’  “  ^\e,Ni^ZC, 


le,A^j  '*cr5*^le,A^2  ’ 

tT 

le,A^2  ‘ 


(62) 


Each  result  in  Eq.  (62)  must  then  be  point-multiplied  with  before  inverting  to  return  to  the 
sequence  domain.  The  final  convolved  result  represented  by  the  matrix  D  is  thus  expressible  as 
the  sum 

r  ^  r  (63) 

^2e,N^T,sa  ®  ^2e,N^T,ss  ®  ^T^2e,N2- 

Figure  5(a)  depicts  the  resulting  image  from  Eq.  (63).  Figure  5(b)  depicts  the  image  result¬ 
ing  from  a  conventional  DFT  multiplication  to  implement  circular  convolution.  Note  the  simi¬ 
larity  between  the  two  images  and  the  effects  of  blurring  caused  by  atmospheric  turbulence  in 
each.  The  average  absolute  difference  between  the  pixels  in  Figures  5(a)  and  (b)  is  0.024  gray 
levels,  which  is  0.0094%  of  the  total  dynamic  range  of  256  gray  levels.  The  maximum  difference 
between  any  two  pixels  is  0.270  gray  levels,  which  is  0.1 1%  of  the  dynamic  range.  The  very 
slight  differences  between  the  images  is  caused  by  the  very  small  outer  values  in  the  PSF  wrap¬ 
ping  back  into  the  image  in  a  slightly  different  manner  for  circular  and  symmetric  convolution. 
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Figure  5.  Degraded  Satellite  Image  Showing  the  Effects  of  Atmospheric 
Turbulence  on  a  1  m  Circular  Aperture  (Negative  Shown  for  Clarity) 


To  summarize  the  procedure  used  in  this  example  and  to  implement  the  symmetric  convo¬ 
lution-multiplication  property  in  general,  the  following  list  of  steps  is  provided: 

(1)  Decompose  the  filter  impulse  response  into  parts  having  support  only  in  a  single 
quadrant  (or  orthant,  if  D  >  2). 

(2)  Decompose  the  parts  existing  in  a  single  quadrant  (or  orthant)  into  antisymmet¬ 
ric  and  symmetric  parts  about  each  axis. 

(3)  Add  parts  with  like  symmetry  from  each  decomposition,  retain  only  the  first- 
quadrant  values,  and  pad  with  zeros. 

(4)  Apply  a  trigonometric  transform  to  each  dimension  of  the  data. 

(5)  Select  appropriate  transform  relations  from  Table  3  based  on  the  type  of  sym¬ 
metry  selected  in  step  (1)  and  the  transforms  selected  in  step  (4). 


59 


(6)  Apply  the  forward  transforms  to  the  impulse  response  matrices  calculated  in 
step  (3). 

(7)  Point-multiply  each  part  of  the  decomposed  and  transformed  impulse  response 
with  the  transformed  data  from  step  (4). 

(8)  Apply  the  inverse  transforms  determined  from  step  (5)  to  each  point-multiplied 
transform-domain  sequence. 

(9)  Add  the  results  of  the  inverse  transforms  to  yield  the  symmetrically  convolved 
sequence. 

In  this  chapter,  the  results  of  deriving  each  of  the  forty  forms  of  the  symmetric  convolu¬ 
tion-multiplication  property  for  discrete  trigonometric  transforms  has  been  presented  by  showing 
how  the  transforms  diagonalize  a  matrix  which  represents  the  symmetric  convolution  operation. 
Derived  in  this  manner,  the  symmetric  convolution-multiplication  property  extends  easily  to 
multiple  dimensions.  The  filtering  of  multidimensional  asymmetric  sequences  is  then  possible 
because  symmetric  convolution  is  equivalent  to  multiplication  in  the  transform  domain  for  each 
of  the  underlying  types  of  symmetry  in  an  asymmetric  image.  The  next  step  in  this  development 
is  to  seek  filters  in  the  trigonometric  transform  domain  which  remove  the  effects  of  blurring 
caused  by  PSFs  which  are  similar  to  the  PSF  used  in  this  example. 


60 


IV.  Image  Reconstruction  Using 
Symmetric  Convolution 

The  theory  developed  in  this  chapter  applies  the  symmetric  convolution-multiplication 
property  of  the  discrete  trigonometric  transforms  (DTTs)  [34]  to  the  linear  image  reconstruction 
problem  [10].  The  presentation  includes  a  derivation  of  one  and  two-dimensional  inverse  and 
scalar  Wiener  filters  expressed  in  the  trigonometric  transform  domain.  For  finite  sequences, 
point-wise  multiplication  in  the  trigonometric  transform  domain  is  equivalent  to  symmetric  con¬ 
volution  in  the  sequence  domain.  Previous  applications  of  the  discrete  cosine  transform  (DCT) 
to  linear  image  reconstruction  [26],  [27]  provided  very  good  diagonal  approximations  for  certain 
types  of  covariance  matrices.  The  DCT  cannot,  however,  simultaneously  diagonalize  the  matrix 
representing  degradation  in  the  linear  model.  Now,  with  the  symmetric  convolution- 
multiplication  property,  very  good  scalar  filter  approximations  are  possible  because  an  exact  di- 
agonalization  of  the  degradation  matrix  is  achievable,  while  still  retaining  the  near  optimum  ap¬ 
proximation  to  the  diagonal  form  of  the  object  covariance  matrix.  In  this  chapter,  the  specific 
forms  needed  from  among  the  forty  cases  of  symmetric  convolution  are  first  reviewed,  and  then 
the  derivations  of  trigonometric  one  and  two-dimensional  inverse  and  scalar  Wiener  filters  are 
presented.  The  next  section  also  presents  a  subtle  property  on  the  equivalence  of  symmetric  con¬ 
volution  to  linear  convolution. 

4.1  Equivalence  Between  Symmetric  Convolution  and  Linear  Convolution 

Throughout  the  remainder  of  this  dissertation,  only  four  of  the  forty  one-dimensional  cases 
of  symmetric  convolution  from  Table  3  are  considered.  This  section  presents  these  cases  along 
with  the  sixteen  two-dimensional  cases  which  result  from  applying  the  four  cases  to  the  rows  and 
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columns  of  matrices.  In  the  latter  two  sections  of  this  chapter,  these  specific  cases  are  used  to 
derive  inverse  and  scalar  Wiener  filters  which  recover  an  object  from  distorted  data.  The 
equivalence  between  symmetric  and  linear  convolution  for  appropriately  zero-padded  sequences 
[34],  [35]  is  first  developed  in  this  section.  This  equivalence  between  symmetric  and  linear  con¬ 
volution  is  similar  to  the  equivalence  which  exists  between  circular  and  linear  convolution  which 
also  arises  from  appropriate  zero-padding  in  the  sequence  domain  [37]. 

The  remaining  theory  presented  in  this  dissertation  requires  only  four  of  the  forty  one¬ 
dimensional  cases  of  symmetric  convolution  from  Table  3.  These  cases,  expressed  in  vector- 
matrix  notation,  are 


(64) 

(65) 

(66) 

and 

dl  =  S2lM[CuA®S2e,N0\ 

(67) 

The  vector  6  represents  a  one-dimensional  sequence  9{n)  which  is  finite  and  zero  outside  the 
interval  0  <  n  <  - 1.  It  has  underlying  half-sample  symmetry  in  both  the  left  and  right  direc¬ 

tions  (HSHS)  in  Eqs.  (64)  and  (66)  and  underlying  half-sample  antisymmetry  in  both  the  left  and 
right  directions  (HAHA)  in  Eqs.  (65)  and  (67).  Recall  that  the  underlying  symmetry  for  a  trigo¬ 
nometric  transform  is  similar  to  the  underlying  periodicity  required  by  the  discrete  Fourier  trans¬ 
form  (DFT).  The  vectors  h[  and  hi  represent  the  right  halves  of  impulse  responses  which  are 
whole-sample  antisymmetric  and  whole-sample  symmetric  in  both  the  left  and  right  directions 
(WAWA  and  WSWS),  respectively.  These  impulse  response  right  halves,  A^and/tJ,  distort  the 
vector  6  to  produce  the  data  sequences  d^,  dl,  d^,  and  </J.  The  subscripts  'a'  and  's'  refer  to 
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antisymmetric  and  symmetric  parts  and  the  prime  superscript  distinguishes  between  data  vectors 
calculated  using  different  cases  of  symmetric  convolution. 

As  previously  demonstrated,  symmetric  convolution  requires  only  the  right  half  of  a  filter's 
impulse  response,  because  the  symmetry  is  required  instead  of  implied.  If  the  sequence  repre¬ 
senting  the  impulse  response  of  the  system  causing  distortion,  h(n),  is  not  symmetric  initially,  it 

must  be  decomposed  into  its  antisymmetric  and  symmetric  parts,  hl{n)  and  hl{n),  which  are  rep¬ 
resented  by  the  vectors  h'  and  AJ.  Then  either  Eq.  (64)  or  (65)  can  convolve  the  symmetric  part 
and  either  Eq.  (66)  or  (67)  can  convolve  the  antisymmetric  part  [34],  [35].  The  four  cases  in 
Eqs.  (64)  -  (67)  are  all  based  on  the  type-II  DCT  for  even  length  sequences,  represented  by  the 
matrix  C2e,N-  The  other  matrix  transforms  in  Eqs.  (64)  -  (67)  are  C,g  ,  and  which 

represent  the  type-I  DCT,  the  type-I  discrete  sine  transform  (DST),  and  the  type-II  DST,  respec¬ 
tively,  all  for  even-length  sequences. 

A  matrix  multiplication  operation  can  perform  symmetric  convolution  in  the  same  way 
that  a  circulant  matrix  performs  circular  convolution.  Multiplications  by  the  diagonal  matrices 

=  diagliSig  and  =  diag|Cig  ;^rAj|  replace  the  point-wise  multiplications  in 

Eqs.  (64)  -  (67).  The  subscript  T' refers  to  the  trigonometric  transform  domain.  These  di- 
agonalizations  result  in  the  expressions 


~  ^2e,N^T,a^2e,N^  ~ 

(68) 

^'a  ~  ~^2e,N^T,a^2e,N^  -  ^SC,a^^ 

(69) 

-  ^2e,N‘^T^^2e,N^  “  ^SC,s^’ 

(70) 

and 

-  ^2e,N^T^^2e,N^  ~  ^'sC,s^- 

(71) 
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The  matrices  H^sc,a>  H^c^s  symmetric  convolution  matrices  for  the  four 

types  of  symmetric  convolution  of  interest  here.  It  is  important  to  note  that  even  though  ^ 

and  Hgc^  *  with  appropriate  zero-padding  in  the  sequences  h^,  h^,  and  9  as  de¬ 
scribed  in  [34]  and  [35],  then  and  d^  =  d^.  The  equality  of  this  vector  representation  is 

assured  because  of  the  equality  which  exists  for  the  sequences  d^(n)  =  d^(n)  and  d^(n)  -  d'^n) 
as  demonstrated  by  Martucci  [34].  Appropriate  zero-padding  assures  that  the  symmetric  convo¬ 
lution  of  the  two  vectors,  h  and  6,  will  equal  the  result  from  linear  convolution.  The  equiva¬ 
lence  of  symmetric  and  linear  convolution  is  completely  analogous  to  the  equivalence  of  circular 
and  linear  convolution  which  also  results  from  appropriate  zero-padding  [37].  Thus,  provided 
the  vectors  and  hi  have  the  correct  underlying  symmetry  and  all  sequences  are  appropriately 
zero-padded,  the  results  of  applying  different  types  of  symmetric  convolution  in  Eqs.  (68)  -  (71) 
are  the  same  because  they  equal  the  result  from  linear  convolution. 

The  equivalence  between  S5anmetric  and  linear  convolution  applies  in  two  dimensions  as 
well.  Consider  the  lexicographically-ordered  vector,  6,  representing  the  Aj  x  Aj  object  matrix, 
0.  Applying  a  type-II  DCT  to  both  dimensions  represented  within  the  vector  9  yields 

®  ^2e,Vj  (^2) 

in  the  transform  domain.  The  subscript  'ss'  shows  that  the  matrix  0  has  underlying  half-sample 
symmetry  about  both  and  «2  in  the  sequence  domain  (HSHS-HSHS).  Equation  (72)  is  the 

lexicographic  equivalent  of  =  C2e^N,^2e,N^  used  for  the  example  in  the  last  section  of  the 
previous  chapter. 

Allow  the  Aj  X  A2  matrix,  H,  which  can  be  viewed  as  a  point  spread  function  (PSF),  to 
represent  the  other  sequence  to  be  convolved.  Recall  that  the  matrices  H^,  H",  HU,  and  HU 
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represent  the  decomposition  of  H  into  its  symmetric  and  antisymmetric  parts.  The  superscript 
'rr'  refers  to  the  PSF's  right  halves  about  both  and  n^.  If  the  vectors  hZ,  C,  C,  and  hZ 

represent  lexicographic  reorderings  of  the  above  decomposed  components  of  then  the  correct 
transformations  to  apply  are 


^T,aa  -  ®  ^le,N^  )Ka  ’  ^r,as  -  ®  ) 


(73) 


based  on  the  different  types  of  symmetry  about  the  and  Wj  axes.  The  expressions  in  Eq.  (73) 
are  lexicographically  equivalent  to  the  expressions  in  Eq.  (62)  because  each  of  the  DTTs  is  sepa¬ 
rable  [27].  Each  result,  4^^^,  4^^^,  and  4^,^,  fromEq.  (73)  must  then  be  point-multiplied 
with  the  vector  .9j.  „  before  inverting  to  return  to  the  sequence  domain.  If  the  lexicographically- 
ordered  vector  d  represents  the  final  convolved  result,  then 


which  is  lexicographically  equivalent  to  Eq.  (63). 

Just  as  in  the  one-dimensional  case,  the  N^N2  x  NiN2  diagonal  matrices 


(74) 


'^T,aa  =  diag{^"^„  },  =  diagj^"^}, 

=diag{4;;„},  and  = 

can  replace  the  point-wise  multiplications  in  Eq.  (74).  The  A^i772  ^  ■^1-^2  symmetric  convolution 
matrices  for  the  two-dimensional  case  are  thus 

^BSC,aa  -  (‘^2e,V,  ®  ^2e,N2  )'^T,aaif-'2e,N,  ®  ^2e,Ni  )’  (^6) 

^BSC.as  -  {f'2e,Ni  ®  ^2e,N2  ^‘^Tfis{f'2e,Ni  ®  ^2e,N2  )’  (^7) 
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and 


^BSC,sa  ~  {^2e,Ni  ®  ^2e,N^  ]'^T,sa  {^2e,N,  ®  ^2e, Af,  )’  (78) 

^BSC^s  -  (^2e,Af,  ®  ^2e.N2  )^7’.ij(^2e,W,  ®  ^2e,JV2  )’  (79) 

where  the  subscript  'BSC  indicates  the  matrices  perform  block  symmetric  convolution  in  two 
dimensions. 

Notice  that  Eqs.  (76)  -  (79)  apply  the  convolution  rules  from  Eqs.  (68)  and  (70),  but  these 
convolution  rules  are  not  unique  since  Eq.  (69)  produces  the  same  result  as  Eq.  (68),  and  Eq.  (71) 
produces  the  same  result  as  Eq.  (70),  as  long  as  appropriate  zero-padding  exists  for  all  sequences. 
Just  as  two  representations  yielded  the  same  result  in  each  of  the  two  one-dimensional  cases, 
there  are  now  four  representations  which  yield  the  same  result  for  each  of  the  four  two- 
dimensional  cases,  producing  a  total  of  sixteen  cases  in  all.  The  transform  relation  acting  on  ei¬ 
ther  the  rows  or  the  columns  in  two  dimensions  may  have  an  equivalent  form  from  Eqs.  (68)  - 
(71)  substituted  for  it.  Thus  the  family  of  two-dimensional  symmetric  convolution  matrices  for 
the  antisymmetric-antisymmetric  portion  of  the  decomposed  PSF  is: 

^BSC,aa  -  (‘^2e.V,  ®  *^2 J.Vj  )^r,aa  {f'2e,N^  ®  ^26,^2  )»  (76) 

^BSCfia  -  ”(^2e.W|  ®  *^26,^2  ^‘^T,aa{^2e,N^  ®  ^2e,N^  )>  (^9) 

mSC,oa  =  -(■Sri.V,  ®C-,lN)‘»T^a[C2e,N,  ®  -S2e,V2).  (81) 

•l^SSC.aa  ~(^2e,Vi  ®  ^26,^2  (*^26,^1  ®  *^26,^2  )•  (82) 

The  family  for  the  antis5anmetric-symmetric  portion  is: 

^BSCfis  -  {f'2l,N^  ®  *^26,^2  ')'^T^{f'2e,Ny  ®  ^2e,N2  )’  (^"^) 

^ BSC, as  ~  (‘^2eVi  ®  *^26,^2  (‘^2e,A^,  ®  ^2e,N2  )»  (8^) 

^BSC,as  -  “(^2e.V,  ®  ^2e,N2  ^‘^Tjas{f'2e,Ni  ®  ^2e,N2  )’  (8^) 
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^BSCflS-  (‘^2e,A',  ®  ^26,^2  )^r.aj(‘^2e,Wi  ®  *^26,^2  )•  (85) 

The  family  for  the  symmetric-antisymmetric  portion  is: 

^^BSC,sa  -  (*^2e,Afi  ®  ^2e,JV2  )‘^T,sa{f^2e,Ni  ®  ^2e,N^ )»  (78) 

^BSC^a  ~  ~{f'2e,N^  ®  ^2e,N^  )^7'^i3  (*^26, A/,  ®  ^2e,N^ )’  (86) 

^BSC,sa  ~  {^2e,Ni  ®  ^2e,N2  )'^T,sa{f'2e,Ni  ®  ^2e,N^ )»  (87) 

^BSC^a  -  “(C'2i,Ar,  ®  ^2e,N^  )'^7^a  (*^26,^,  ®  ^2e,N^  )■  (88) 

The  family  for  the  symmetric-symmetric  portion  is; 

^BSC,ss  -  {^2l,Ni  ®  ^2l,N2  ')‘^T^s{f'2e,N^  ®  ^2e,N^  )j  (79) 

^BSC^s  -  (‘^2e,Wi  ®  ^2e,W2  )^Oi(‘^2e,W,  ®  ^2e,N^  )»  (89) 

^BSC.ss  -  {^2e.N,  ®  ^2e,N2  )'^T,ss{^2e,N,  ®  ^2e,N^  )’  (90) 

^BSC.ss  ~  {^2e,Ni  ®  *^26,^2  )^7',M(‘^2e.Ar,  ®  ^2e,N^  )•  (91) 

Here  again  none  of  the  symmetric  convolution  matrices  within  a  particular  family  are  exactly 

equal,  i.e.  ^  ^SSC^a  ^  ^BSC,aa  ^BSC,aa’  ^ BSC, as  ^  ^ BSC, as  ^  ^BSCfis  ^ BSC, as’ 

^ BSC, sa  ^  ^BSC,sa  ^  ^BSC,sa  ^BSC,sa  ’  ^BSC,ss  ^  ^BSC,ss  ^  ^BSC,ss  ^BSC,ss-  However,  with 

appropriate  zero-padding,  the  vectors  =  Hgscfla^,  d'^a=H'Bsc^a^’  =HBscjaa^’  and 
d'aa  =  H's'sc^a^  will  all  be  equal;  the  vectors  =  HBsc,as^,  d;^  =  H’Bsc,as^,  C  =  H'^sc,as^’  and 
das=HBsc,as^  will  all  be  equal;  the  vectors  d,„  =  Egsc^a^,  d'^  =H'Bsc^„e,  rf"  =  H'^sc,sa^’  and 
C  =  ^BSc.sa^  will  all  be  equal;  and  the  vectors  d,,  =  Hbsc,ss^,  d',,  =  H'gsc^s^,  d';,  =  H'^sc,ss^’ 
and  d"J = H'bsc  ss^  will  all  be  equal.  The  resulting  sequences  will  also  equal  the  sequence  which 
arises  from  circular  convolution  with  appropriate  zero  padding,  because  all  are  equal  to  the  result 
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from  linear  convolution.  The  same  vector  d  will  result  if  any  of  the  different  equivalent  forms  in 
Eqs.  (76)  -  (91)  are  substituted  into  Eq.  (74),  as  long  the  vectors  h  and  6  are  appropriately  zero 
padded.  All  of  the  different  forms  of  symmetric  convolution  in  Eqs.  (76)  -  (91)  use  the  vector- 
matrix  form  [9]  of  Martucci's  [34]  symmetric  convolution-multiplication  property  for  trigono¬ 
metric  transforms  derived  in  the  previous  chapter. 

The  equivalence  of  symmetric  convolution  to  linear  convolution  as  outlined  in  this  section 
holds  in  the  trigonometric  transform  domain  for  any  linear  shift-invariant  filtering  application. 
These  results  are  applied  in  the  following  sections  by  recasting  some  traditional  linear  image  re¬ 
construction  filtering  operations  into  the  trigonometric  transform  domain. 

4.2  Inverse  Filtering  in  the  Trigonometric  Transform  Domain 

The  results  of  the  previous  sections  are  valid  for  any  finite  sequences  represented  by  the 
vectors  d,  h,  and  0,  which  are  all  lexicographically-ordered  for  the  two-dimensional  case.  The 
only  underlying  assumption  was  appropriate  zero-padding  to  demonstrate  the  equivalence  of  the 
convolutional  forms.  In  this  section  the  vector  d  represents  a  detected  image  which  arises  from 
blurring  an  original  object  vector  ^  by  a  point  spread  function  (PSF)  represented  lexicographi¬ 
cally  as  h.  The  goal  of  the  inverse  filters  derived  in  this  section  is  to  recover  a  vector  estimate  of 
the  object,  6,  from  the  vector  </ in  the  trigonometric  transform  domain.  This  classical  problem 
as  it  applies  to  image  reconstruction  finds  the  two-dimensional  impulse  response  of  a  filter,  rep¬ 
resented  lexicographically  by  the  vector  f  which  recovers  0  from  d,  given  knowledge  of  h.  This 
is  the  same  problem  which  generated  Fourier  transform  domain  solutions  presented  in  Sec¬ 
tion  2.2  as  background.  The  following  subsections  present  equivalent  representations  in  the 
trigonometric  transform  domain,  first  for  one  and  then  for  two  dimensions. 
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4.2.1  The  One-Dimensional  Inverse  Filter  for  Trigonometric  Transforms.  The  devel¬ 


opment  of  this  subsection  derives  a  result  in  the  trigonometric  transform  domain  similar  to 
Eq.  (10)  which  presented  the  one-dimensional  inverse  filter  in  the  Fourier  domain.  Using  the 

same  notation  as  before  whereby  d  =  H9  and  6  =  Fd  =  FHd,  understand  that  in  this  case  the 
matrices  H  and  F  represent  one-dimensional  symmetric  convolution  matrices  rather  than  circular 
convolution  matrices.  The  notation  here  drops  the  subscript  'SC  which  was  useful  in  previous 
sections  to  distinguish  between  symmetric  and  circular  or  skew-circular  convolution.  The  dis¬ 
cussion  from  here  forward  concerns  itself  exclusively  with  symmetric  convolution  and  trigono¬ 
metric  transforms.  The  'BSC  subscripts  will  likewise  not  appear  in  two-dimensional  block  sym¬ 
metric  convolution  matrices,  nor  will  the  subscript  T  appear  in  trigonometric  transform  domain 
quantities. 

The  symmetric  convolution  matrices  Ef  and  Fmay,  in  general,  be  asymmetric,  which  re¬ 
quires  a  decomposition  into  their  symmetric  and  antisymmetric  parts  so  that  d  =  {H^  +H^6  and 

e  =  {F,  +  F,)d  =  (F,  +  F,){H,  +  H,)e.  (92) 


Expanding  Eq.  (92)  and  making  substitutions  from  Eqs.  (68)  -  (71)  results  in 


^  —  |(“^2e,v5a‘^2e.Af  )(‘^2e,V^a^2c.V  )  "'■  (  ^le.N^a^le.N  )(‘^2e,V^Ae.V  ) 

{f'2e,N^s^2e,N  )(“^2e,V^o*^2e,V  )  "*■  {f'2e,N^s^2e,N  )(^2e,V^i^2e,Af  )]^- 


(93) 


Equation  (93)  incorporates  substitutions  of  the  different  equivalent  forms  of  symmetric  convolu¬ 
tion  based  on  the  assumption  that  the  sequence  domain  data  vector,  d,  is  the  result  of  linear  con¬ 
volution.  This  assumption  implies  that  sufficient  zero-padding  [34],  [35]  exists  in  the  vectors  h 
and  6  to  ensure  their  equivalence.  The  innermost  matrices  inside  each  term  of  Eq.  (93)  multiply 

to  identity,  so  that  exact  reconstruction  to  yield  the  vector  0  =  9,  requires 

7s'»s-’Pa'»a=h  (94) 
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and 


(95) 


where  /  is  an  identity  matrix,  and  0  is  an  Ny.'N  zero  matrix.  Because  all  the  matrices  in 

Eqs.  (94)  and  (95)  are  diagonal,  the  equations  are  equivalent  to  the  matrix  equation 

for  k  =  0,  1,  ...  ,  A^-1.  The  terms  1?(,(^),  9a{k),  and?^(A:)  inEq.  (96)  represent  the 

A:-Ath  terms,  of  the  matrices  and^^,  respectively. 

Solving  Eq.  (96)  for  and  P^k)  yields 


(96) 


and 


9M  = 


^l{k)  +  P‘^{k) 


(97) 


^sik) 

P‘l{k)  +  P‘^{ky 


(98) 


Equations  (97)  and  (98)  are  the  one-dimensional  inverse  filter  expressed  in  the  trigonometric 
transform  domain  for  symmetrically  convolved  sequences.  If  the  sequence  h(n)  possesses 
strictly  whole-sample  symmetry,  i.e.,  h^in)  =  then  the  inverse  filter  reduces  to 


7s{k)  = 


1 

■ 


(99) 


Equation  (99)  is  similar  to  the  Fourier  domain  result  in  Eq.  (10).  The  following  subsection 
shows  similar  results  for  two  dimensions. 

4.2.2  The  Two-Dimensional  Inverse  Filter  for  Trigonometric  Transforms.  In  two  di¬ 
mensions,  H  and  F  become  two-dimensional  block  symmetric  convolution  matrices  which  may 
also  in  general  be  asymmetric.  They  therefore  require  a  decomposition  into  their  symmetric  and 
antisymmetric  parts  about  both  the  and  «2  axes  so  that  d  =  and 
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(100) 


^  -  {^aa  + 

=  {^aa  +  Fas  +  Psa  +  ^ss){Haa  +  +  H^^)d. 

Expanding  Eq.  (100)  and  making  substitutions  from  the  equivalent  forms  in  Eqs.  (76)  -  (91)  re¬ 
sults  in 


®Qiv, 


®  ^le,N^  )[(*^2e,JV,  ®‘52.U)  ®C^e.N^) 

+  (‘^2e.W|  ®‘52tvJ^?'a.(‘^2..V.  ®C2e.N,) 

■*■  (‘^2e,V,  ®  *^2e,A^2  ]^sa[^2e,Ni  ®  ^2e,N2  ) 

®‘^U)  ^ss{^2e,Ni  ®‘5'2e.V,)] 
®  ‘y2..vJ[(C  2e,Ni  ®  S;iN)‘^aa[S  2e,N^ 

~{f'2e,N^  ®  ^2l,N^‘^as{f'2e,N^  ®^2e,N^ 
+  (^2e,V|  ®>y2‘eU)  ^sa{^^2e,N^ 

-(QU  ®  ®  -52.,^,)] 

®  ^2e.W2)[(‘^2e,V,  ®  ^2e,N^‘^aa{^2e,N^  ®  ^2e,N^ 

®C2U)  ^as{^2e,N^  ®^2e,N2) 

~{^2e,Ni  ®^2e,N2)'^sa{f'2e,Ni  ®^2e,//2) 
“(‘^2e,A^i  ^ss{^2e,Ni 

®C2e,V2)[(C  2e,Ny  ®  ^2e,N2)'^aa[^2e,N^  ®  *^2e,A^2) 
”  {^2e,N^  ®  ^2e,N2  )'^as{f'2eMi  ®  ^2e,N2  ) 
“  (^eVi  ®  ^2e,N2  )'^sa  {^2eMi  ®  ^2e,A^2  ) 
^ss{f'2e,Ny 


e. 


(101) 


Equation  (101)  incorporates  substitutions  of  the  different  equivalent  forms  of  the  convolution 
matrices  in  Eqs.  (76)  -  (91).  The  substitutions  are  permissible  because  of  the  assumption  that  the 
lexicographically-ordered  data  vector,  d,  arises  as  the  result  of  linear  convolution  in  the  sequence 
domain. 
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For  this  two-dimensional  case,  the  interior  matrices  in  each  term  of  Eq.  (101)  again  multi¬ 


ply  to  identity,  so  that  exact  reconstruction  requires 


"Paa^aa  ~  Pastas  ~  PsQp^sa  +  Ps^ss  -  ^ (102) 

PaaP^as  +  PasP^aa  ~  Psa'^ss  ~  PssP^sa  =  0»  (103) 

PaaP^sa  ~  PasP^ss  +  PsaP^aa  ~  PssP^as  =  0,  (1 04) 

and  PaaP‘ss+PasP^sa+PsaP‘as+PA=^-  (105) 


Each  matrix  in  Eqs.  (102)  -  (105)  is  diagonal,  which  generates  the  matrix  equation 


~P‘aa{KA)  -P^asiK^K)  I^^XKA)' 

'PaaAA) 

T 

P‘asiKK)  P^aaiKh)  -P^saAA) 

Pas  A  A) 

0 

P^saAA)  ~P^ssAA)  P^aaAA)  ~P‘asAA) 

PsaAA) 

0 

P^ssAA)  P‘saAA)  P^asAA)  P‘aaAA)  . 

PssAA). 

_0_ 

for  ki  =0,  1,  ...  ,  N^-l  and^Tj  =  0,  1,  ...  ,  N2-I.  The  terms  P^asi^u^i)’ 

P^sa(.K^k2),  PaaikM,  Pasi^M,  and  ?„(^„/:2)  in  Eq.  (106)  repre¬ 
sent  the  diagonal  elements  of  the  7/jA^2^^i-^2  niatrices  Pas-’  Psa-> 


and  respectively.  Equation  (106)  has  the  solution 


PaaikM  = 


,  ^  )  I  -P^aa(kiA  )P‘li.K  A)  +  '^P^asiKA  )P‘saiKh  )P‘ssiK ,  ^  )J 


(107) 


Pasi^l>^)- 


-1 

^(^lA) 


[ P^L  (^1 5  ^  )  +  P^as  (^1  >  ^  )P^m  (^1 5  ^  )  “  P^as  (^1 » ^2  )‘^sa  (^1  >  ^2  )  | 

1  +  P^asi.KA)P^l(.Kk2)  +  '^P‘aa{KA)P^sa(.kA)P^ss(kM] 


,  (108) 


and 


PsaiK’K)  = 


5ji(^i>^)  ~ 


-1  P‘L(K’^)'^P‘sa(^l’^2)P‘aa(^l’^)  P‘sa(^’^)P‘asi^\’^)  1 

■^(^15^)  I  +  P‘sa(K’^)P^ss(.^\’^)'^^P‘aaih’^)P^asih’^)P‘ss(^\’^)\ 

1  \p^ssi^l>^)~P^ss(^l’^)P^aai^l’^)'^P^ssi^l’^)P^as(^l’^)  1 


,  (109) 


(110) 


where  the  determinant  of  the  4x4  matrix  in  Eq.  (106)  is  given  by 
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(Ill) 


Equations  (107)  -  (1 1 1)  are  the  two-dimensional  inverse  filter  for  trigonometric  transforms.  If 
the  sequence  /z(«„«2)  possesses  strictly  whole-sample  symmetry  about  the  and  «2  axes,  i.e., 

Kainuih)  =  ^aa(ki,k2)  =  0,  A„,(«i,«2)  =  =  0,  and  //,„(«!, ^2)  =  ^sa(ki,k2)  =  0,  then 

the  two-dimensional  inverse  filter  reduces  to 


9ss(K,k2) 


1 


(112) 


Equation  (112)  closely  resembles  its  Fourier  domain  equivalent  in  Eq.  (12). 

The  trigonometric  transform  domain  realizations  of  the  inverse  filter  for  symmetrically- 
convolved  one  and  two-dimensional  sequences  expressed  in  Eqs.  (97),  (98),  and  (107)  -  (1 1 1) 
suffer  from  the  same  high-frequency  gain  problem  which  plagued  their  Fourier  domain  equiva¬ 
lent  forms  for  noisy  sequences.  The  problem  becomes  quite  serious  when  the  image  model  ex¬ 
pands  to  incorporate  noise  with  a  uniform  power  spectral  density  across  all  frequencies.  In  the 
next  section,  a  method  is  presented  of  regularizing  the  high  frequency  gain  of  the  inverse  filter  in 
the  presence  of  noise. 


4.3  Wiener  Filtering  in  the  Trigonometric  Transform  Domain 

In  this  section,  a  derivation  of  the  scalar  Wiener  filter  is  presented  in  the  trigonometric 
transform  domain  [12].  A  Wiener  filter  introduces  a  degree  of  regularization  to  the  inverse 
problem  and  is  more  capable  of  filtering  data  models  with  noise,  as  explained  in  the  background 
section  on  Wiener  filtering  in  Chapter  II.  To  incorporate  noise  into  the  data  model  of  the  previ¬ 
ous  section,  the  model  must  have  an  added  term  so  that  it  now  becomes  d  ~  HO  +  w.  The  matrix 
jfiT  is  a  symmetric  convolution  matrix  for  the  one-dimensional  case  and  a  block  symmetric  convo- 
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lution  matrix  for  the  two-dimensional  case.  The  vector  w  is  again  a  zero-mean  uniform-variance 
noise  vector  whose  samples  are  uncorrelated  both  with  the  samples  of  the  object  vector,  0,  and 
with  each  other.  The  vectors  d  and  9  represent  finite  sequences  for  one  dimension  and  are  lexi¬ 
cographically-ordered  representations  of  finite  matrices  for  two  dimensions.  The  object  vector, 
9,  is  itself  random  with  constant  mean  vector  Hq  =  9,  which  can  be  assumed  to  equal  0  without 

loss  of  generality  [28],  and  correlation  matrix  ^dd  =  99^ .  The  Wiener  filter  seeks  the  vector 

yv  -  ^ 

estimate,  9,  of  9  which  minimizes  the  mean  squared  error,  e,  where  the  vector  s=9-9. 

The  solution  to  this  problem  from  [28]  appeared  previously  for  circulant  and  block  circu- 
lant  matrices  in  the  Fourier  transform  case.  The  equivalent  to  Eq.  (14)  for  symmetric  convolu¬ 
tion  is 


9  = 

which  produces  the  recovery  filter 

f^r,,h'^\hr,,hUr,J[\ 


(113) 


(114) 


Recall  that  the  background  section  on  Wiener  filtering  in  Chapter  II  described  vector  and  scalar 
Wiener  filters  in  one  and  two  dimensions  expressed  in  the  Fourier  transform  domain.  In  that 
case  the  matrices  F  and  H  were  circulant  for  one-dimensional  filters  and  block  circulant  for  two 
dimensional  filters.  Here  in  the  trigonometric  case,  the  matrices  F  and  ^  will  represent  sym¬ 
metric  convolution  matrices  for  one-dimensional  filters  and  block  symmetric  convolution  matri¬ 
ces  for  two-dimensional  filters. 

The  Fourier-domain  scalar  Wiener  filters  presented  for  background  in  Eqs.  (23)  and  (27) 
for  one  and  two  dimensions  under  the  assumption  of  wide-sense  stationarity  for  the  object,  pro¬ 
vide  a  good  approximation  to  the  vector  or  generalized  Wiener  filter  [38].  The  approximation 
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arose  because  the  discrete  Fourier  transform  did  not  exactly  diagonalize  the  symmetric  Toeplitz 
form  of  the  correlation  matrix  for  a  wide-sense  stationary  object  [48].  Note  that  the  noise  is  al¬ 
ready  stationary  because  it  has  zero  mean  and  its  samples  are  uncorrelated  with  each  other. 

In  the  derivations  which  follow  in  the  trigonometric  transform  domain,  the  data  vector,  d, 
must  arise  as  the  result  of  a  linear  convolution  of  the  object  vector,  0,  with  the  point  spread 
function,  h.  This  assumption  of  linear  convolution  allows  the  substitution  of  different  equivalent 
forms  of  the  two-dimensional  convolution  matrices  from  Eqs.  (76)  -  (91)  into  the  components 
which  result  from  decomposing  the  matrix  H  in  Eq.  (114).  In  the  one-dimensional  problem,  the 
one-dimensional  convolution  matrices  from  Eqs.  (68)  -  (71)  may  be  substituted.  This  assumption 
of  linear  convolution  also  exists  in  Fourier-domain  derivations  because  the  processes  which  gen¬ 
erate  blurred  data  are  linear  and  not  circular  in  nature.  Circular  convolution  arises  because  it  is 
the  underlying  form  of  convolution  for  DFTs  and  it  is  convenient  mathematically  to  process  the 
data  using  DFTs.  However,  for  the  data  to  arise  from  linear  convolution,  the  object  must  be 
zero-padded,  but  a  zero-padded  or  support-constrained  sequence  cannot  be  wide-sense  stationary. 
The  assumption  of  linear  convolution  underlying  the  process  which  generated  the  blurred  data 
thus  seems  to  contradict  the  assumption  of  wide-sense  stationarity  for  the  object. 

Hunt  and  Cannon  [25]  addressed  this  conundrum,  and  their  work  is  further  refined  in  [51]. 
Their  work  chooses  the  more  accurate  nonstationary  object  model  with  support  constraints  to 
account  for  the  zeros  which  must  exist  to  be  equivalent  to  linear  convolution.  A  model  without 
support  constraints  does  not  reflect  the  true  statistics  of  the  object.  The  correlation  matrix  of  an 
object  with  support  constraints  has  zeros  which  appear  due  to  the  zero-padding  which  must  exist 
for  linear  convolution  to  appear.  These  zeros  do  not  appear  for  an  object  without  support  con¬ 
straints.  A  model  without  support  constraints  is,  however,  more  mathematically  tractable  be¬ 
cause  it  is  well-approximated  by  a  diagonal  matrix  in  the  transform  domain  [48].  In  a  support- 
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constrained  model,  a  higher  degree  of  correlation  exists  between  frequencies  in  the  Fourier  do¬ 
main  [36].  The  purpose  of  this  research  is  to  find  scalar  filters  which  require  good  diagonal  ap¬ 
proximations  in  the  trigonometric  transform  domain.  In  this  case  it  is  better  not  to  use  support 
constraints  in  the  model  so  that  good  diagonal  approximations  result.  For  vector  filters,  the  sup¬ 
port-constrained  nonstationary  object  model  is  better  because  it  exhibits  enhanced  performance 
by  incorporating  the  off-diagonal  correlations  into  the  filter  design  [14].  The  choice  of  a  model 
without  support  constraints  is  also  more  appropriate  in  this  scalar  filter  derivation  because  this 
effort  is  the  first  attempt  to  apply  the  symmetric  convolution-multiplication  property  to  a  Wiener 
filter.  The  original  Fourier-domain  Wiener  filters  were  all  unconstrained  [15],  [27],  [28],  [43], 
[48],  so  this  choice  provides  a  better  comparison  to  earlier  work.  In  the  following  discussion,  the 
one-dimensional  scalar  Wiener  filter  for  trigonometric  transforms  is  derived  first,  and  then  its 
two-dimensional  equivalent  is  presented. 

4.3.1  The  One-Dimensional  Scalar  Wiener  Filter  for  Trigonometric  Transforms.  The 
first  step  in  deriving  the  one-dimensional  scalar  Wiener  filter  is  to  substitute  the  decomposed 
symmetric  convolution  matrices  F  =  Fg+F^  and  H  =  Hg+H^  into  Eq.  (114).  Making  these 
substitutions  and  then  bringing  the  bracketed  expression  over  to  the  left  side  produces 

+f;)[(^„  ^h,)r^{h,  +h,Y + rJ^ = R^{H,  (115) 

This  derivation  then  follows  a  similar  procedure  as  the  inverse  filter  derivation.  It  expands 
Eq.  (115)  and  then  substitutes  the  equivalent  forms  from  Eqs.  (68)  -  (71)  to  produce 


{f^2e.N^a^2e,N\^2l,N^a^2e,N)^ee{^2e,N^a^2e,N) 
'^{f^2e,N^a^2e,N  ^^2e,N^J-2e.N  {f'2e,N^s^2e,N  ) 

'^{^^2e,N^a^2e,N  \^2e,N^s^2e,N)^O0{^2e,N^a^2e,N') 

'^{^2e,N^a^2e,N  )[^2l,N^s^2e,N  )^0a{^2l,N^s^2e,N  ) 
'^{^2e,N^s^2e,N  )(‘^2e,V^a^2e,v)-^w(*^2i.V^a^2e,w)  •  •  • 
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+(‘^2eV?i‘^2e,w)(‘^2e,Ar'^aQe,Af)-®ee(^e,W^iQe,w) 

+(  ^e,N^s  Qe,  W  )(  ^e,N^s^2e,N)^0e  {^2e,N^a  ^2e.  W  ) 
+(^i,w5^i^2e,Ar)(^2e,W^i^2e.A^)-^6'(  ^2e,N^s^2e,N) 
'^{~^2l,N^a^2e,N)^w^L,N^2lN  +  {^l,N^s^2e,N)^ww^2e,N^2l 

~  ^2l,N^2e,N^09{^2e,N'^a^2e,N)  +  ^2e,N^2e,N^e0{f'2e,N'^s^2e,N^ 


-T 

■'2e,N 


Equation  (116)  simplifies  to 


)  +  9s^s^0,0,'^a\^2l. 
+C2lM[-9a^aR0^0^^s  +  9s('»sR0,0,9^s  + 
+‘^2e.Jv[5a^s^0,0,^a  '^9s‘^a^0,0,‘^a 


-T 

2e,N 


S~^ 

^2e,N 

-f-r 

■'2e,N 


+S2lN[9a9^s^0^0^'»s  +  9s‘^aR0,0:»s 
~  ^2e,N^0,0j^a^2e,N  ^2e,N^0^9‘s^2e,N  ’ 


(116) 


(117) 


where  —  C2e,N^e0^2e,N’  -  ^2e,N^ww^2e,N’  ^^Tn,  -  ^2e,N^ww^2ejf  ■  The  Simpli¬ 
fication  of  Eq.  (1 16)  to  produce  Eq.  (1 17)  uses  the  facts  that  =  f?(f  and  .  For 


Eq.  (117)  to  hold,  it  must  follow  that 


-9a{9‘o^0,0,9^a  ^W^7ff,)  +  9s9^s^0,0,9‘a  -R©,0,^a. 

(118) 

-9a"^a^0,0  9^s  +  9s{^st^0,0  9^s  +  =  ^0fi^s> 

(119) 

9a9^sR0^0:»a  +  9s9fa«~0^0^a  =  «. 

(120) 

and 

9.9fsRe^e:»s  +  9s9fa^0^0:»s  =  «■ 

(121) 

Every  matrix  in  Eqs.  (1 18)  -  (121)  is  either  exactly  or  approximately  diagonal.  Recognize 
that  the  matrices  “Pa^  9 s,9f a,  and  are  always  diagonal  because  they  are  the  diagonalized 

forms  of  symmetric  convolution  matrices  in  the  trigonometric  transform  domain.  The  matrices 
and  R^^  will  be  diagonal  because  the  samples  of  w  are  assumed  to  be  uncorrelated  with 
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each  other  and  to  have  uniform  variance,  ,  so  that  =  a^l.  In  general,  the  transform- 
domain  correlation  matrix  will  be  well-approximated  by  its  diagonal  elements  if  the  sam¬ 
ples  of  the  object,  6,  are  highly  correlated  with  each  other  [27].  Specifically  if  Rgg  is  the 
N'kN  correlation  matrix  for  a  Markov-I  process,  [38],  [47],  [53],  then  the  correlation  matrix 
in  the  DCT  domain  will  approximate  a  diagonal  matrix.  For  a  Markov-I  process,  the 

m  -  «th  element  of  R^  is  =  ycJ'"  "I  where  Q)<  p<\.  As  an  example  of  how  well  the 

transform  domain  correlation  matrix  approximates  a  diagonal  matrix,  let  Rgg  be  the  256  x  256 
correlation  matrix  of  a  Markov-I  process  with  p  =  0.9.  Then  the  matrix  Rq^@^  =  n 

in  the  DCT  domain  has  98.75%  of  its  total  energy,  i.e.,  the  sum  of  the  squares  of  its  elements, 
along  the  diagonal.  As  a  comparison  the  Fourier  domain  matrix  =  W^^RggW^^  has 

98.20%  of  its  total  energy  contained  in  its  diagonal  elements.  The  transform  domain  correlation 
matrix,  R@^0^  ,  will  approach  an  exactly  diagonal  matrix  either  as  N  increases  without  bound  or 

as  p  increases  to  1 .  This  approximate  diagonalization  occurs  because  the  DCT  provides  an  ex¬ 
cellent  approximation  to  the  eigenvectors  of  the  correlation  matrix  of  a  Markov-I  process  [27], 
[39].  Thus  under  the  above  assumptions  for  the  noise  and  object  processes,  R^^  andJL,^.  are 

exactly  represented  by  their  diagonal  elements  70^ {k)  and  Tofik)  for  k  =  0,  1,  ...  ,N,  and 
is  well-approximated  by  its  diagonal  elements  6>J(^)  for  ^  =  0,  1,  ...  ,  N -\. 

The  above  approximations  allow  for  a  scalar  solution  based  solely  on  the  diagonal  ele¬ 
ments  of  the  matrices  in  Eqs.  (1 18)  -  (121).  Solving  for  the  diagonal  elements  of  Eqs.  (118)- 
(121)  produces  the  4x2  overdetermined  matrix  equation 
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€^s(k) 

-l^aikWk) 

9a(k) 

-9‘ai.k) 

&]{k) 

9sik) 

^a(kWk) 

^a(k) 

0 

^a(k)1^s(k) 

0 

(122) 


Equation  (122)  reduces  to 


7  7  TV 

0l 


0 

0 


0?  0? 


(123) 


where  each  term  in  Eq.  (123)  depends  on  the  index,  k,  but  temporarily  has  the  indexing  sup¬ 
pressed  to  preserve  space  in  the  equation.  The  system  of  equations  will  be  consistent,  and 
Eq.  (123)  will  equal 


9a(k) 

-9^a(k) 

0]{k) 

9s(k) 

9^sik) 

0]{k)_ 

(124) 
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0,  ^  =  0 

7(/Xk)  =  -2N(T\  k  =  \,2,  (126) 

4iVcr^  k  =  N 


4Na^,  k  =  0 

2N(j\  k  =  \,2,  N-l  (127) 

0,  k  =  N 

whenever  =  cr^I.  From  Eqs.  (126)  and  (127)  it  follows  that  Eq.  (125)  holds  for 

^  =  1,  2,  It  is  also  straightforward  to  verify  Eq.  (125)  at  A:  =  0  by  using  the  fact  that 

=  0  based  on  its  required  symmetry  in  the  transform  domain  [34].  Verifying  that 


Eq.  (125)  holds  at  k  =  N  requires  the  facts  that  O^iN)  =  0  and  =  0  based  on  their  re¬ 

quired  symmetiy  in  the  transform  domain  [34].  Then  L'Hopital's  rule  must  be  applied  to  take  the 
limit  as  @1{N)  ^  0  of  the  derivative  of  the  numerator  and  denominator  of  each  side  of  Eq.  (125) 
with  respect  to  0\{N). 

Thus  Eq.  (125)  is  valid  for  all  values  of  k,  and  Eq.  (124)  will  hold,  from  which  matrix 
multiplication  followed  by  scalar  division  produces 

9a(.k)  = - ,  (128) 

^ 
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and 


9sik)^ 


(129) 


'^sik) 


0]{k) 


Equations  (128)  and  (129)  are  the  one-dimensional  scalar  Wiener  filter  expressed  in  the  trigono¬ 
metric  transform  domain  for  symmetrically  convolved  sequences.  Note  first  that  for  the  noise- 


free  case  where  w  =  0,  which  implies  that  T&Xk)  =  0  and  Ti/^k)  =  0,  Eqs.  (128)  and  (129)  re¬ 
duce  to  Eqs.  (97)  and  (98),  and  the  scalar  Wiener  filter  becomes  an  inverse  filter.  Note  also  that 
if  the  sequence  h(n)  possesses  strictly  whole-sample  symmetry,  i.e.,  h„{n)  =  =  0,  then  the 

scalar  Wiener  filter  reduces  to 


m)= 


^sik) 


0l{k) 


(130) 


Equation  (130)  closely  resembles  Eq.  (23)  for  the  one-dimensional  Fourier  case.  Similar  results 
for  two  dimensions  are  shown  in  the  following  subsection. 

4.3.2  The  Two-Dimensional  Scalar  Wiener  Filter  for  Trigonometric  Transforms.  The 
derivation  of  the  two-dimensional  scalar  Wiener  filter  in  the  trigonometric  transform  domain  is 
somewhat  more  involved  than  it  was  for  one  dimension.  The  two-dimensional  version  of  the 
general  asymmetric  cases  for  the  filter  matrices  requires  a  four-way  decomposition  for  each 
block  symmetric  convolution  matrix,  so  that  F  =  F^„+F^  +  F^^+  F^  and  H  =  + 

Making  these  substitutions  into  Eq.  (1 14)  yields 


{Ka  +Ks+  ^sa  +  +  ^as  +  +  ^ss)^9e{^aa  +^as+  ^sa  +  ^  ^ 

^  ^ee{^aa+^as+^aa+^as)  ■ 

The  expansion  of  Eq.  (131)  produces  68  terms  on  the  left-hand  side.  The  Kronecker  product 
transformations  of  the  convolution  matrices  in  Eqs.  (76)  -  (91)  must  then  transform  each  of  these 
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terms  into  the  trigonometric  transform  domain.  The  two-dimensional  versions  of  Eqs.  (116)- 
(121)  will  subsequently  become  quite  cumbersome.  The  details  of  this  lengthy  derivation  appear 
in  the  appendix.  The  end  result  in  the  appendix  is  a  16  x  4  overdetermined  system  of  equations 
similar  to  Eq.  (122),  which  reduces  to  the  7x4  matrix  equation 


^2  "^aa 

^aa  ^ 

-^aa'^as 

'»aa'»ss 

"Paa 

P‘aa 

^as  ^  :r 

0^ 

^ss 

^as^sa 

-^as'^ss 

“Pas 

-P‘as 

-^aa^sa 

'^as'^sa 

^sa  :r 

-"^sa^ss 

Psa 

^aa^ss 

^as  ^ss 

^sa  ^ss 

1  Tip' 

Pss 

^as 

^aa 

-^sa 

0 

-^ss 

^aa 

-^as 

0 

^ss 

^as 

^aa 

0 

- 

All  the  quantities  in  the  matrix  equation  are  indexed  over  A:,  and  ytj,  but  Eq.  (132)  temporarily 
suppresses  the  indexing  to  conserve  space. 

Unfortunately  the  set  of  equations  in  Eq.  (132)  is  not  consistent  as  it  was  in  the  one¬ 
dimensional  case.  Thus  no  general  solution  exists  in  the  trigonometric  transform  domain  for  the 
general  case  of  a  PSF,  A(n,,n2),  which  possesses  all  four  types  of  underlying  symmetry.  Solu¬ 
tions  exist  for  /(«,,n2)  in  the  sequence  domain  from  solving  Eq.  (131)  directly,  but  this  involves 
a  large  matrix  inversion  for  even  small  image  and  PSF  sizes.  Transform  domain  representations 
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are  computationally  easier  to  implement,  so  the  following  discussion  focuses  on  the  conditions 
under  which  a  solution  to  Eq.  (132)  exists. 


The  solution  7aai^’K)-  =  0  exists  trivially  for  the 

case  when  =  =  0,  but  is  of  little  use.  A  solution 

also  exists  for  the  four  cases  where  any  three  of  the  four  components  are  zero.  If  ‘»^^(ki,k2)  = 

~  ~  ^  ^00(^15^2)^®*  then  the  PSF  A(«|,«2)  would  possess  strictly 

whole-sample  antisymmetry  about  both  the  and  rij  axes,  and  the  scalar  Wiener  filter  for  this 
case  is 


9aai.Kk2)  =  - 


(133) 


®ssik\,k2) 

For  the  case  where  =  0  and  ^  0,  the  scalar 


Wiener  filter  is 


9asiki,k2) 


‘^as(kljk2) 


'^as(ki,k2)  + 


'^as(.k\,k2) 

0l(kuk2) 


(134) 


For  the  case  where  '»^iki,k2)  =  ??^„,(A:,,^2)  =  “^ssikM  =  0 


and  7‘^a(.k\,k2)  ^  0,  the  scalar 


Wiener  filter  is 


9saikx,k2)  =  - 


-9‘sa{K,k2) 


^ssik\,k2) 


(135) 


For  the  case  where  '»^{k^,k2)  =  7f„^{k^,k2)  =  '»,„(^i,^)  =  0  and  7f,,{k\,k2')  ^  0,  the  scalar 


Wiener  filter  is 
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(136) 


&]s{kM 


Equation  (136)  is  the  scalar  Wiener  filter  for  the  case  of  a  PSF  which  possesses  strictly  whole 
sample  symmetry  about  both  the  n■^  and  axes. 


Note  that  the  subscript  '55'remains  in  the  term  @]Xk^,k2)  for  all  four  cases  in  Eqs.  (133)  - 
(136).  Throughout  the  derivation  of  the  two-dimensional  scalar  Wiener  filter  presented  in  the 
appendix,  substitutions  of  equivalent  convolution  matrices  in  the  transform  domain  from 
Eqs.  (76)  -  (91)  has  resulted  in  a  transform  domain  object  correlation  matrix  of  the  form 

^0,,0ss  ~  (^2e.w,  ®  ^2e,N2)^eo{f^2e,Ni  ®  ^2e,Afj )  •  Choosing  transform  relations  so  that  the  trans¬ 
form  domain  object  correlation  matrix  always  lies  in  the  transform  domain  of  the  type-II  DCT  is 
what  allows  it  to  be  well-approximated  by  its  diagonal  elements  and  yield  a  scalar  filter. 

Equation  (132)  will  also  produce  a  consistent  set  of  equations  in  four  of  the  six  cases 
where  two  of  the  four  types  of  underlying  symmetry  are  simultaneously  nonzero.  The  equations 
will  be  inconsistent  for  the  two  cases  where: 

(0  "^aa (^1 » ^2  )  “  ^ssik\  >  ^2  )  ~  ‘^asi.k\ » ^2  )  ^  ^sa (^I » ^2  )  ^ 

and  (//)  ^ai(^i5^2)  ~  ^00(^15^2)  ^  and 

They  will  be  consistent  for  the  four  cases  where: 

(0  ^aa(^l>^2)  ~  ^fli(^l’^2)  ~  (^1 » ^2 )  ^  ^^^(^1 5^2 )  ^ 

(«)  =  0,  ^  0>  ^nd  0; 

(“0  =  '^ssiki,k2)  =  0,  '»aaiki,k2)  ^  0,  and  ^?^^<,(A:,,^2)  ^  0; 

and  (iv)  (k^ , ^2 )  =  (^i > ^2 )  =  0,  l^aa (^1 , ^2 )  and  ,k2)^0. 
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For  the  first  case  where  =  =  and 

the  resulting  scalar  Wiener  filter  will  be 


(137) 

9as(kl,k2)  =  0, 

(138) 

9sa(Jh,k2)  = - ^ ^ 

9‘l{Kk2)  +  9‘l{k„k2)  + 

&]si.K,k2) 

(139) 

and  7ss{k„k2)  = - ^iKA) - 

&li.K,k2) 

(140) 

For  the  second  case  where  '»,AkM  =  I^.AkM  =  0,  ^^^kM  ^  0,  and  ^.XkM  ^  0,  the 

resulting  scalar  Wiener  filter  will  be 

‘Paaikx,k2)  =  (}, 

(141) 

9as(kM  = - Z^^sikM  ^ ^ 

9^i(kM  +  1i^l(k„k2)+ 

^ssikM 

(142) 

9sa(ki,k2)  =  0, 

(143) 

and  9ss(k„k2)  = - ^ - 

9‘iik„k2)  +  1i‘l{k„k2)  +  ^^^^ 

&l(kM 

(144) 

For  the  third  case  where  ‘^^XkM  =  9‘^(.K,k2)  =  0,  l^aaikM  ^  0,  and  ^^^aiki, 

^2)  ^  0,  the 

resulting  scalar  Wiener  filter  will  be 

9‘l(.K,k2)  +  1^l{k„k2)  +  ^^'^^^ 

Oli-kM 

(145) 
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9as(kl,k2)  =  0, 


(146) 


9sai.KA)  = - ^ 

'^laikM  +  KikM  +  ^^s^^ 

^ssiKA) 

and  ?,,(^„^2)  =  0.  (148) 

For  the  fourth  case  where  "^.aikM  = ‘»ss(kM  =  0,  ‘»„a(k^,k2)*0,  and  '»^Xki,k2)  ^  0,  the 
resulting  scalar  Wiener  filter  will  be 

%a(kl,k2)  = - ^aaiklA)  ^ 

^aa(kM  +  ^l(k„k2)  +  ^^^ 

^ssikM 

%sikM  = - Z:^aXkM  ^ 

0l(kM 

9sa(ki,k2)  =  0,  (151) 

and  %Xkvk2)  =  0.  (152) 

Note  that  for  each  of  the  eight  nontrivial  cases  where  the  two-dimensional  scalar  Wiener 
filter  exists,  it  reduces  to  the  inverse  filter  of  Eqs.  (107)  -  (1 1 1)  whenever  the  noise  terms  go  to 
zero.  The  set  of  equations  in  Eq.  (132)  is  inconsistent  for  each  of  the  four  cases  having  three 
nonzero  terms  and  for  the  general  case  of  four  nonzero  terms. 

Two  options  are  available  to  find  a  solution  to  the  general  problem  of  a  PSF  having  all 
four  types  of  underlying  symmetry.  The  first  option  finds  the  least-squares  solution  to  the  over¬ 
determined  system  in  Eq.  (132).  The  second  option  makes  different  substitutions  from  equiva¬ 
lent  forms  for  the  transforms  to  yield  a  4  x  4  system  of  equations. 
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Using  shorthand  notation  to  express  the  7x4  matrix  equation  in  Eq.  (132)  as  “,¥1  =  the 
least-squares  solution  is  ^  The  inverse  of  exists  because  has  rank 

four.  The  vector  is  the  solution  which  minimizes  ||^  -  [46],  Unfortunately,  both  the 

least  squares  solution  and  the  solution  to  the  4  x  4  system  of  equations  which  results  from 
choosing  different  equivalent  forms  for  the  transforms  are  too  lengthy  to  be  presented. 

Each  of  these  two  methods  of  finding  a  general  solution  has  other  disadvantages  as  well. 
The  disadvantage  to  finding  the  least  squares  solution  to  Eq.  (132)  is  that  it  is  not  exact  and  will 
therefore  introduce  more  error  in  addition  to  the  errors  caused  by  noise  and  by  the  scalar  ap¬ 
proximation  which  retains  just  the  diagonal  elements  of  the  object  correlation  matrix.  The  dis¬ 
advantage  to  choosing  different  transforms  in  the  derivation  of  the  two-dimensional  filter  is  that 

the  result  depends  on  the  terms  6>^„(A:„^2).  and  which  are  all  not  as 

well-approximated  by  their  diagonal  elements  as  eilXKh)  [26].  The  energy  in  the  off-diagonal 

terms  which  is  not  accounted  for  in  the  diagonal  approximation  will  introduce  additional  error  as 
well. 

The  fact  that  no  exact  general  solution  exists  to  Eq.  (132)  is  not  as  restrictive  as  it  might 
first  appear.  The  form  of  the  solution  appearing  in  Eq.  (136)  based  on  a  PSF  having  only  whole- 
sample  symmetry  is  the  form  of  the  filter  which  will  likely  prove  the  most  useful  in  practice,  be¬ 
cause  whole-sample  symmetry  is  a  necessary  condition  for  a  filter  to  have  linear  phase  [37]. 

Many  image  reconstruction  methods  model  asymmetric  PSFs  with  inherent  nonlinear  phase,  such 
as  those  used  to  model  the  effects  of  atmospheric  turbulence  [40],  as  random  processes.  In  these 
cases,  even  though  individual  realizations  of  a  PSF  are  asymmetric,  the  mean  PSF  will  often  be 
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symmetric  so  that  only  whole-sample  symmetry  will  be  present  in  the  scalar  Wiener  filter. 
Equation  (136)  is  a  valid  form  of  the  scalar  Wiener  filter  in  these  cases. 

The  existence  of  inverse  and  scalar  Wiener  filters  in  the  trigonometric  transform  domain 
for  one  and  two  dimensions  has  been  demonstrated  in  this  chapter.  The  two  dimensional  scalar 
Wiener  filter  is  limited  by  the  type  of  symmetry  present  in  the  PSF  in  certain  cases.  It  was 
claimed  that  the  two-dimensional  whole-sample  symmetric  version  is  the  most  useful  case  for 
image  processing.  In  the  following  chapter,  the  performance  of  the  scalar  Wiener  filter  for  this 
particular  case  is  analyzed. 
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V.  Performance  of  the  Scalar  Wiener  Filter 
for  Trigonometric  Transforms 

The  performance  of  the  two-dimensional  scalar  Wiener  filter  is  analyzed  in  this  chapter 
under  the  assumption  that  the  degradation  system  point  spread  function  (PSF)  is  whole-sample 
symmetric.  Systems  possessing  this  characteristic  are  the  most  likely  to  be  encountered  in  prac¬ 
tice.  Systems  with  random  asymmetric  PSFs  often  have  a  mean  PSF  which  is  whole-sample 
symmetric.  To  analyze  the  scalar  Wiener  filter's  performance,  an  example  is  first  provided  and 
then  the  normalized  mean-squared  error  is  calculated  for  several  PSFs  and  objects  with  varying 
parameters.  The  chapter  concludes  with  a  brief  discussion  on  the  computational  advantages  of 
the  trigonometric  transform-based  filter. 

5.1  An  Example 

In  this  section,  the  performance  of  the  new  trigonometric  transform  versions  of  the  two- 
dimensional  inverse  and  scalar  Wiener  filters  is  demonstrated  by  providing  an  example.  It  is 
shown  how  a  blurred  object  can  be  completely  recovered  with  an  inverse  filter.  An  inverse  filter 
cannot,  however,  recover  the  original  object  in  the  presence  of  noise.  A  scalar  Wiener  filter  then 
recovers  an  estimate  of  the  object  from  its  noisy  blurred  version  by  regularizing  the  high- 
frequency  gain  of  an  inverse  filter.  The  results  of  this  example  compare  an  estimate  of  the  object 
using  a  trigonometric  scalar  Wiener  filter  to  an  estimate  using  a  traditional  Fourier  scalar  Wiener 
filter.  This  example  serves  to  visually  demonstrate  the  effects  of  image  reconstruction  filters  in 
the  trigonometric  transform  domain. 

The  256  x  256  pixel  satellite  object  shown  in  Figure  4  for  the  trigonometric  transform 
domain  filtering  example  in  Section  3.4  forms  the  basis  of  this  example  as  well.  Figure  6  shows 
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the  effects  of  inverse  filtering  in  the  trigonometric  transform  domain  for  a  blurred  noiseless  ob¬ 
ject  and  a  blurred  noisy  object.  In  Figure  6(a),  the  object  is  blurred  using  a  16x16  pixel  Gaus- 


c.  Noisy  Blurred  Object  (SNR  =  20  dB)  d.  Noisy  Object  Unrecoverable  with 

Inverse  Filter 


Figure  6.  Inverse  Trigonometric  Filtering  Example  Using  Hanning- Windowed 
Gaussian  Degradation  Filter  with  Me  Width  of  4  Pixels 
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sian-shaped  PSF  windowed  with  a  Hanning  window  [37]  to  prevent  an  abrupt  cutoff  at  the  tran¬ 
sition  point.  The  shape  of  the  PSF  rolls  off  to  a  value  of  1/e  at  a  distance  of  4  pixels  from  the 
center.  Applying  the  two-dimensional  inverse  filter  of  Eq.  (1 12)  in  the  trigonometric  transform 
domain  results  in  the  image  shown  in  Figure  6(b).  The  original  object  is  completely  recovered 
from  its  blurred  version  based  on  knowledge  of  the  blurring  PSF.  Figure  6(c)  shows  the  results 
of  adding  noise  with  a  signal-to-noise  ratio  (SNR)  of  20  dB  to  the  blurred  object.  After  adding 
noise,  the  object  is  no  longer  recoverable  with  an  inverse  filter  as  demonstrated  by  the  result  in 
Figure  6(d).  The  reason  the  trigonometric  inverse  filter  no  longer  recovers  the  object  from  the 
noise  is  that,  like  the  Fourier  inverse  filter,  it  amplifies  high-frequency  components  of  the  noise. 

Figure  7  presents  the  results  of  applying  Fourier  and  trigonometric  scalar  Wiener  filters  to 
the  noisy  blurred  object  of  Figure  6(c).  Equation  (27)  implements  the  scalar  Wiener  filter  in  the 


Figure  7.  Scalar  Wiener  Filtering  Example  Using  Same 
Degradation  Filter  as  Inverse  Filtering  Example 
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Fourier  transform  domain,  and  Eq.  (136)  implements  the  scalar  Wiener  filter  in  the  trigonometric 
transform  domain.  Unlike  the  inverse  filters  for  Figure  6,  the  scalar  Wiener  filters  for  Figure  7 
use  a  block  processing  technique  for  their  implementation.  The  need  for  a  block  technique  arises 
because  the  correlation  matrices  which  regularize  the  inverse  filter  become  very  large,  even 
though  a  scalar  Wiener  filter  is  only  concerned  with  their  diagonal  elements.  These  diagonal 

elements  in  the  Fourier  domain  are  the  and  terms  in  Eq.  (27).  The  diago¬ 

nal  elements  in  the  trigonometric  transform  domain  are  the  terms  7(^l{k^,k2)  and  ©lXkx,k2)  in 
Eq.  (136).  The  block  processing  technique  divides  the  noisy  blurred  image  of  Figure  6(c)  into 
16x16  blocks,  zero  pads  each  block  to  a  size  of  32  x  32,  performs  all  calculations  in  32  x  32 
tiles  of  transform  space,  inverse  transforms  the  results  back  to  the  spatial  domain,  and  then  sums 
the  overlapping  results.  Since  a  Wiener  filter  depends  on  an  object  to  originate  from  a  random 
process  with  known  mean  and  covariance,  this  example  uses  the  statistics  of  the  object  itself  to 
calculate  the  matrices  and  Cgg.  The  object  has  a  nonzero  mean  which  yields  a  covariance 
rather  than  a  correlation  matrix.  The  noise  variance  is  set  equal  to  1/1 00th  of  the  object  variance 
to  achieve  a  20  dB  SNR. 

Comparing  Figures  7(a)  and  (b),  it  is  clear  that  the  trigonometric  scalar  Wiener  filter  pro¬ 
duces  a  better  quality  estimate  of  the  object  than  the  Fourier  scalar  Wiener  filter.  One  criticism 
of  Wiener  filtering  is  that  it  tends  to  oversmooth  an  image  to  compensate  for  noise.  This  effect  is 
clearly  more  pronounced  in  the  Fourier  case,  because  the  object  estimate  is  still  quite  blurry  al¬ 
though  the  noise  has  been  averaged  out.  The  trigonometric  case  achieves  a  much  sharper  con¬ 
trast  and  higher  spatial  frequency  resolution  because  of  the  improved  energy  compaction  of  the 
DCT.  Comparing  the  results  of  Figures  7(a)  and  (b)  to  the  original  object,  a  mean-squared  error 
normalized  to  the  total  number  of  pixels  is  345.5  for  the  Fourier  case  and  21 1 .7  for  the  trigono- 
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metric  case.  Expressed  as  a  ratio  to  the  square  of  the  highest  gray  scale  value  of  256,  these  be¬ 
come  -22.7  dB  and  -24.9  dB,  respectively.  These  values  fall  within  the  ranges  of  some  more 
general  expressions  for  the  mean-squared  error  of  the  Fourier  and  trigonometric  scalar  Wiener 
filters,  as  will  be  seen  from  the  results  of  the  following  section. 

5.2  Mean-Squared  Error  Performance 

Closed-form  expressions  exist  for  the  mean-squared  error  of  Wiener  filters  because  Wie¬ 
ner  filtering  is  a  linear  technique.  Using  the  method  of  [28]  and  [38],  the  mean-squared  error  is 

expressed  as  Ijflj  =  ||,  which  expands  to 

\sf^=lr[Ree-2FHRee+F[HR,eHUR,J)F^].  (153) 

The  transform  domain  expression  for  Eq.  (153)  is 

||rfg  =  Tr{i?^0 -29^R0e  +  (154) 

where  the  matrix  T represents  either  a  DFT  or  a  DCT  matrix.  All  other  quantities  in  Eqs.  (153) 
and  (154)  have  been  defined  previously.  It  is  easier  to  calculate  the  mean-squared  error  using 
Eq,  (154)  in  the  transform  domain  because  all  of  the  matrices  are  diagonal  except  the  transform 
domain  correlation  matrix,  Rqq.  Recall  that  the  diagonal  transform  domain  filter  matrix,  “P,  de¬ 
pends  on  a  diagonal  approximation  of  Rqq.  The  terms  in  the  diagonal  approximation  are  the 

0\{k^,k2)  terms  from  the  Fourier  domain  case  of  Eq.  (27),  and  the  €P'^fk^,k2)  terms  from 
Eq.  (136)  for  the  trigonometric  case.  To  accurately  assess  the  error  that  this  approximation  in¬ 
troduces,  the  error  expression  must  compare  the  diagonalized  version  imbedded  in  the  filter  ma¬ 
trix,  P,  to  the  nondiagonal  version  in  Eq.  (154)  which  it  approximates.  Each  transform  has  dif¬ 
ferent  scale  factors  which  cause  different  gains  in  the  transform  domain.  The  error  calculations 


93 


in  the  transform  domain  must  account  for  these  different  scale  factors.  The  error  calculations 


must  then  apply  the  appropriate  scale  factors  to  calculate  the  mean-squared  error  which  results 
back  in  the  spatial  domain  for  both  the  Fourier  and  the  trigonometric  transforms. 

Figures  8-11  display  the  results  of  calculating  the  error  using  Eq.  (154)  in  various  two- 
dimensional  filtering  scenarios.  Each  filtering  scenario  depicted  uses  Hanning-windowed  Gaus- 
sian-shaped  PSF  filters  to  model  the  degradation  similar  to  the  previous  example.  The  scenarios 
adjust  parameters  for  the  length  (A^  of  the  Nx  N  PSF,  the  correlation  coefficient  (p)  in  the 
Markov-I  object  correlation  matrix,  the  signal-to-noise  ratio  (SNR),  and  the  width  in  pixels  to  the 
point  at  which  the  PSF  has  a  value  of  l/e.  Adjusting  these  parameters  creates  different  scalar 
Wiener  filters  in  Eqs.  (27)  and  (136)  for  the  Fourier  and  trigonometric  transform  domains.  Each 
graph  displays  the  mean-squared  error  (MSE)  for  the  Fourier  scalar  Wiener  filter  with  a  solid 
line,  and  the  MSE  for  the  trigonometric  scalar  Wiener  filter  with  a  dashed  line. 

In  Figure  8,  the  PSF  filter  dimension,  N,  increases  while  the  correlation  coefficient,  p,  is 
set  to  0.9.  Additionally  the  SNR  is  fixed  at  20  dB,  and  the  PSF  \/e  width  increases  at  a  constant 


Figure  8.  Normalized  Mean-Squared  Error 


vs.  Filter  Dimension 
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rate  of  Ml  6  pixels  as  increases.  In  each  of  the  remaining  Figures  9  - 1 1,  the  filter  dimension, 
N,  is  held  to  16x16  pixels.  Figure  9  shows  how  the  MSB  decreases  as  p  increases  for  PSF  1/e 
widths  of  1  and  8  pixels,  with  the  SNR  fixed  at  20  dB.  Figure  10  shows  the  MSB  decreasing 
with  increasing  SNR  for  values  of  p  of  0.5,  0.9,  and  0.99,  with  the  PSF  1/e  width  fixed  at  1  pixel. 


Figure  9.  Normalized  Mean-Squared  Brror 
vs.  Correlation  Coefficient 


Figure  10.  Normalized  Mean-Squared  Brror 
vs.  Signal-to-Noise  Ratio  (SNR) 


Finally  Figure  1 1  shows  the  MSB  increasing  as  the  filter  1/e  width  increases  and  the  filter  band¬ 
width  decreases.  The  figure  shows  curves  for  SNRs  of  0, 20,  and  40  dB,  with  p  fixed  at  0.9. 
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Figure  11.  Normalized  Mean-Squared  Error 
vs.  Point  Spread  Function  (PSF)  Me  Width 


In  all  of  the  cases  tested,  the  trigonometric  scalar  Wiener  filter  performed  better  than  the 
Fourier  scalar  Wiener  filter,  often  demonstrating  an  improvement  in  mean-squared  error  on  the 
order  of  20  -  35%.  The  greatest  increase  in  performance  occurs  as  the  bandwidth  of  the  degrada¬ 
tion  filter  in  the  linear  model  decreases.  The  larger  gap  between  the  DFT  and  DCT  curves  at  the 
higher  SNRs  in  Figure  1 1  demonstrate  this  effect.  The  fact  that  the  gap  is  greater  in  Figure  9  for 
the  case  where  the  PSF  1/e  width  is  8  pixels  than  it  is  for  the  case  where  it  is  1  pixel  also  dem¬ 
onstrates  this  improvement.  The  increased  performance  results  from  the  fact  that  a  type-II  DCT 
possesses  a  near  optimum  energy  compaction  property  about  the  low-frequency  indices  for 
highly  correlated  images  [27].  The  improvement  with  increasing  correlation  appears  in  Figures  9 
and  10. 


5.3  Computational  Complexity 

In  addition  to  its  improved  MSB  performance,  the  trigonometric  scalar  Wiener  filter  has 
the  advantage  of  requiring  fewer  calculations  to  implement.  The  advantage  to  performing  filter¬ 
ing  in  the  trigonometric  transform  domain  is  that  for  real  sequences,  the  transform  domain  coef¬ 
ficients  are  also  all  real,  where  they  are  complex  in  the  Fourier  domain.  Fast  algorithms  exist  for 
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the  DCT  based  entirely  on  real  arithmetic  which  operate  with  the  same  number  of  floating  point 
operations  as  fast  Fourier  transforms  [5],  [27],  [39],  These  real-arithmetic  algorithms  for  the 
DCT  amount  to  a  tremendous  savings  in  computational  performance  since  complex  additions 
require  twice  the  number  of  floating  point  operations  as  real  additions  and  complex  multiplica¬ 
tions  typically  require  six  times  the  number  of  floating  point  operations  as  real  multiplications. 
Fast  hardware  realizations  of  these  trigonometric  image  processing  techniques  are  also  possible 
because  of  their  computational  cost  savings. 

In  this  chapter  the  performance  of  the  scalar  Wiener  filter  has  been  demonstrated  visually 
through  the  use  of  an  example,  and  algebraically  through  calculations  of  the  filter's  mean  squared 
error  performance.  An  argument  has  been  presented  that  the  filter  also  possesses  the  additional 
advantage  of  requiring  fewer  calculations  to  implement.  Thus  the  filter  designed  during  the 
course  of  this  dissertation  research  not  only  has  better  performance  than  its  Fourier  equivalent, 
but  it  is  easier  to  compute  as  well. 
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VI.  Conclusions 


A  summary  of  the  major  results  and  contributions  of  the  research  conducted  for  this  disser¬ 
tation  is  presented  in  this  chapter.  Some  directions  for  future  research  in  this  area  of  study  are 
also  given. 

6.1  Results  and  Contributions 

The  results  of  this  research  demonstrate  how  to  apply  the  recently-developed  symmetric 
convolution-multiplication  property  of  the  discrete  trigonometric  transforms  [34]  to  the  tradi¬ 
tional  image  reconstruction  problems  of  inverse  and  Wiener  filtering.  Each  of  the  forty  forms  of 
the  symmetric  convolution-multiplication  property  for  discrete  trigonometric  transforms  is  shown 
to  exist  as  a  vector-matrix  operation.  The  convolutional  forms  of  the  trigonometric  transforms 
can  then  diagonalize  a  matrix  which  represents  the  symmetric  convolution  operation  [9],  [1 1]. 
Derived  in  this  manner,  the  symmetric  convolution-multiplication  property  extends  easily  to 
multidimensional  asymmetric  sequences  which  represent  the  most  general  type  of  sequences  en¬ 
countered  in  practice. 

The  new  diagonal  forms  for  symmetric  convolution  matrices  represent  the  first  time  such 
forms  have  been  developed  for  the  entire  family  of  all  forty  cases  of  symmetric  convolution. 
Although  diagonal  forms  have  been  shown  to  exist  for  some  of  the  forty  cases  [41],  [42],  and  the 
initial  discovery  of  a  symmetric  convolution-multiplication  property  for  discrete  trigonometric 
transforms  [34]  was  not  part  of  this  research,  the  diagonal  forms  derived  during  the  course  of  this 
research  still  make  a  significant  contribution  to  the  depth  of  knowledge  surrounding  symmetric 
convolution.  As  a  result  of  this  research,  the  relationship  between  the  diagonal  properties  of  dis- 
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Crete  Fourier  transform  matrices  and  discrete  trigonometric  transform  matrices  has  been  related 
for  the  first  time. 

The  diagonalizing  forms  developed  through  this  research  permit  a  clearer  understanding  of 
the  principle  of  how  S3mmetric  convolution  in  the  sequence  domain  is  equivalent  to  multiplica¬ 
tion  in  the  transform  domain.  With  the  new  matrix  formalism  of  symmetric  convolution,  it  is 
much  easier  to  visualize  how  the  symmetric  convolution-multiplication  property  extends  to  se¬ 
quences  in  multiple  dimensions  and  which  might  be  asymmetric.  The  filtering  of  multidimen¬ 
sional  asymmetric  sequences  is  then  possible  because  of  the  equivalence  of  symmetric  convolu¬ 
tion  to  multiplication  in  the  transform  domain  for  each  of  the  underlying  types  of  symmetry  in  an 
asymmetric  image. 

The  research  presented  here  uses  the  newly-derived  vector-matrix  form  of  S5mmetric  con¬ 
volution  to  calculate  for  the  first  time  inverse  and  scalar  Wiener  filters  in  the  trigonometric  trans¬ 
form  domain  [10].  The  new  forms  of  the  inverse  and  scalar  Wiener  filters  closely  resemble  their 
traditional  Fourier  domain  counterparts.  Specifically,  the  new  forms  of  scalar  Wiener  filters  [12] 
are  possible  only  because  of  the  newly-developed  symmetric  convolution  multiplication  property 
of  discrete  trigonometric  transforms  [34].  Previous  applications  of  the  discrete  cosine  transform 
to  linear  filtering  [26],  [27]  provided  very  good  diagonal  approximations  for  certain  types  of  co- 
variance  matrices.  The  trigonometric  transforms  could  not,  however,  simultaneously  diagonalize 
a  matrix  representing  degradation  in  the  linear  model.  Now,  with  the  symmetric  convolution- 
multiplication  property,  very  good  scalar  filter  approximations  are  possible  because  an  exact  di- 
agonalization  of  the  degradation  matrix  is  achievable,  while  still  retaining  the  near  optimum  ap¬ 
proximation  to  the  diagonal  form  of  the  object  covariance  matrix. 

The  two-dimensional  scalar  Wiener  filter  is,  however,  limited  by  the  type  of  symmetry 
present  in  the  point  spread  function  in  certain  cases.  In  general,  a  two-dimensional  sequence  can 


99 


possess  four  types  of  underlying  symmetry.  It  can  have  antisymmetric  and  symmetric  compo¬ 
nents  about  each  of  its  two  axes.  The  two-dimensional  scalar  Wiener  filter  derived  during  the 
course  of  this  research  is  found  to  exist  in  the  trigonometric  transform  domain  only  in  cases 
where  at  most  two  of  the  four  types  of  underlying  symmetry  are  present  at  one  time.  This  limita¬ 
tion  is  not  severe  because  the  two-dimensional  version  of  the  scalar  Wiener  filter  which  has 
whole-sample  symmetry  about  both  axes  is  claimed  to  be  the  most  useful  case  for  image  process¬ 
ing.  The  whole-sample  symmetric  case  is  thus  the  focus  for  analyzing  the  performance  of  the 
scalar  Wiener  filter. 

It  was  demonstrated  that  a  scalar  Wiener  filter  provides  better  mean-squared  error  per¬ 
formance  for  symmetric  point  spread  functions  while  reducing  the  required  number  of  computa¬ 
tions.  The  trigonometric  transform  domain  realizations  of  scalar  Wiener  filters  use  fewer  calcu¬ 
lations  than  traditional  Fourier  realizations  of  the  same  types  of  filters  because  they  use  entirely 
real  arithmetic  for  real  inputs.  The  fact  that  these  filters  provide  better  performance  with  fewer 
calculations  should  prove  beneficial  to  a  wide  variety  of  other  communications,  signal  process¬ 
ing,  and  control  system  applications,  even  though  they  were  derived  specifically  for  image  re¬ 
construction. 

The  need  to  improve  the  quality  of  noisy,  blurred  images  of  spacebome  objects  is  what 
initially  prompted  the  Air  Force  to  sponsor  this  design  of  two-dimensional  trigonometric  scalar 
Wiener  filters.  The  large  number  of  computations  required  for  nonlinear,  iterative  image  recon¬ 
struction  techniques  launched  this  reinvestigation  of  ways  to  improve  traditional  linear  ap¬ 
proaches.  The  results  of  this  research  show  that  the  performance  of  traditional  linear  image  re¬ 
construction  methods  can  be  enhanced  while  also  reducing  the  number  of  computations  required. 
It  is  hoped  that  even  better  performance  can  be  achieved  by  applying  the  symmetric  convolution- 
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multiplication  property  of  discrete  trigonometric  transforms  to  some  of  the  other  more  powerful 
nonlinear  image  reconstruction  techniques  by  reducing  their  number  of  calculations  as  well. 

6.2  Directions  for  Future  Research 

Any  future  research  which  applies  the  symmetric  convolution-multiplication  property  of 
discrete  trigonometric  transforms  to  image  reconstruction  should  attempt  to  compensate  for  the 
assumptions  which  were  required  to  conduct  this  research.  This  research  required  many  as¬ 
sumptions  because  it  was  the  first  time  symmetric  convolution  had  ever  been  used  in  image  re¬ 
construction.  Two  assumptions  which  can  be  eliminated  through  future  research  are  the  assump¬ 
tion  of  an  object-independent  noise  model  and  the  assumption  that  the  point  spread  function 
which  distorts  the  object  is  completely  known. 

The  imaging  model  for  this  research  requires  noise  that  is  zero-mean,  additive,  of  uniform 
variance,  uncorrelated  with  itself,  and  independent  of  the  object.  For  real-world  imaging  situa¬ 
tions,  this  noise  model  is  not  as  accurate  as  one  which  assumes  that  the  noise  is  an  object- 
dependent  Poisson  random  process  [17].  The  filtering  techniques  derived  here  can  be  improved 
by  using  a  more  accurate  Poisson  noise  model.  Incorporating  object-dependent  photon  noise  into 
the  imaging  scenario  would  make  the  techniques  developed  here  compatible  with  more  modem 
image  reconstmction  methods. 

The  second  assumption  which  can  be  eliminated  through  future  research  is  that  the  trigo¬ 
nometric  inverse  and  Wiener  filters  derived  here  require  knowledge  of  the  point  spread  function 
of  the  degrading  system.  In  practice  the  point  spread  function  of  the  degrading  system  is  often 
not  known  exactly.  This  lack  of  knowledge  leads  to  a  technique  called  blind  deconvolution  [4], 
[45].  In  blind  deconvolution  an  estimate  is  sought  for  both  the  object  and  the  unknown  point 
spread  function.  In  blind  deconvolution,  the  point  spread  function  is  unknown  but  deterministic. 
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Two  iterative  techniques  use  inverse  [2]  and  scalar  Wiener  filters  [8]  in  the  Fourier  do¬ 
main  to  estimate  the  object  and  the  unknown  point  spread  function.  Both  of  these  techniques 
should  benefit  from  the  increased  performance  of  trigonometric  inverse  and  scalar  Wiener  filters. 
The  estimates  should  be  of  better  quality  and  the  iterations  should  converge  to  a  solution  faster. 
The  number  of  computations  should  also  be  reduced  because  the  trigonometric  filters  use  real 
arithmetic.  If  future  research  reveals  success  with  these  iterative  techniques,  then  the  symmetric 
convolution-multiplication  property  could  be  applied  to  other  more  advanced  techniques  such  as 
bispectrum  processing  for  multiframe  blind  deconvolution  [40]. 

A  final  area  for  future  research  which  allows  the  point  spread  function  to  be  unknown  is  to 
model  the  point  spread  function  as  a  random  process.  Reasonable  statistical  approximations  to  a 
wide  variety  of  atmospheric  degradations  can  often  be  found  [40].  These  statistical  approxima¬ 
tions  lead  to  wide  variety  of  powerful  techniques  [13],  [14]  which  require  knowledge  of  the  sta¬ 
tistics  of  the  distorting  process  but  not  the  point  spread  function  itself.  These  techniques  all  use 
Fourier  domain  representations  which  could  benefit  from  being  recast  into  the  trigonometric 
transform  domain. 

Theoretically  it  can  be  concluded  from  the  promising  results  of  this  research  that  almost 
any  Fourier-domain  image  reconstruction  technique  should  be  able  to  benefit  from  being  recast 
into  the  trigonometric  transform  domain. 
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Appendix:  Derivation  of  the  Two-Dimensional  Scalar 
Wiener  Filter  for  Trigonometric  Transforms 


Many  of  the  significant  steps  in  the  derivation  of  the  two-dimensional  scalar  Wiener  filter 
for  trigonometric  transforms  are  presented  in  this  appendix.  The  end  result  is  the  7x4  matrix 
equation  displayed  in  Eq.  (132)  of  Section  4.3.2.  The  derivation  starts  with  Eq.  (131)  which  is 
repeated  below 


{Faa  ^Fas  +  F^  +Has+H^  +  + 


(131) 


Recall  that  Eq.  (13 1)  is  the  solution  resulting  from  the  Bayesian  Gauss-Markov  Theorem  [28]  for 
two-dimensional  asymmetric  filters  represented  by  the  matrices  F  =  +  F^^  and 


H  =  -I-  .  The  matrix  H  performs  symmetric  convolution  to  implement  the 

point  spread  function  (PSF)  of  a  system  which  degrades  an  object  represented  by  the  vector  6. 
The  matrix  F  performs  symmetric  convolution  to  recover  a  vector  estimate,  9,  of  the  original 
object. 


Expanding  the  terms  on  the  left-hand  side  of  Eq.  (13 1)  produces 


sa-^^66^^  as  ~  ^  aa^-^  sa-^^09^^  sa 


+FaaHssFeeHL+F,,H,,R,gHl  +  F,,H,^R,,Hi  +F,,H,,R^HI 

'^Fg^HggRggHg^  -f  FgfiggRggHgg  +  Fg^Hg^RggH^g  +  Fg^H gg  RggH 

+Fg,Hg,RggHl  +Fg,Hg,RggHl+Fg,Hg,RggHi-vFgfI^RggHl 
■^FgJHsgReeHlg  +  Fgfl^gRggH^.  +  Fg^H^gRggHjg  +  Fg,H^gRggHj^ 
'^FgsH^gRggHgg  +  Fg^H^^RggHg^  ^  Fg^H^^RggH^g  +  Fg^H^^RggH^^ 

+  FsgHggRggHgg  +  F^g  H gg  RggH g^  +  F^gH gg  RggH ^g  +F^gHggRggHl 

+FgHg^RggHl  +  FM..R»,HL  +  FM..R,ML  +  FM.M..H 


sa'^^  as  '  ^  sa-^-^  as^^OO^"^  sa 


^  sa^^  as'^^OO^'^  ss 
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sa^ee^ aa  ^sa^sa^OO^as  ^sa^sa^dG^sa  "*■  ^sa^sa^GG^ss 
~^^sa^ss^GG^aa  ^sa^ss^GG^as  "*■  ^sa^ss^GG^la  ^  ^sa^ ss^GG^L 
~^^ss^aa^GG^aa  ^ss^ aa^GG^ as ^ss^ aa^GG^ sa  aa^GG^ls 

'^^sJ^as^GG^aa  ^  ^ss^as^GG^as  as^GG^ sa  ^ss^^as^GG^L 

■^Fs,H,,RggHI  +  +  F,,H,^R,,Hl  +  F,fl,^R,,Hl 

'^^aa^ww  ^as^ww  ^sa^ww  ^ss^ww 

=  ^ee^la  +  ^ee^ls  +  ^ee^L  +  ^ee^L  ■  (155) 

After  substituting  the  trigonometric  transform  domain  versions  of  the  symmetric  convolution 
matrices  from  Eqs.  (76)  -  (91),  Eq.  (155)  becomes 

[{C2U  ®<^2U)94^2e.N,  ®C2e,vJ] 

+[(C2.U  ®C2lN,)9aa[S2e,N,  ® ‘S2..V,)][(-S2tv.  ® 

® ‘^2;UK(C’2e,V.  ^<^2..;.,)]" 

+[(C'2j,Af,  ®(^2lM^9aa[S2e,N^  ®  *S'2e,Ar2  )][(‘5'2'i.V,  ®  ^2l,N^'^„^{C2e^^® 

•^^[(‘S2-*V.  ®C-lN,)9‘,a[C2e.,,  ®C2e.M,)f 

+[(<^2eU  ®^2lN,)9aa{S2e,N,  ® ‘y2e,vJ][(‘5'2i,V,  ®  ^2lN,)'^aa{C2e,N,®  C2e,N,)] 

■Ree[(C2lN,  ®C^l2,,)»ss[C2e,N,  ®C2e,N,)\ 

-[(*52;U  ®C2U)9c.[C2e,N,  ®  *52..^  JK^.U  ® ‘^2;'^,  ® 

•^4fc  ®S2U]‘»aa[C2e,N,  ^^2..^,)]" 

“^(‘^2e,yv,  ®^2e,N2)9aa{^2e,Ni  ®  *5'2e,V2)J|(^2e,^/,  ®  ^2e,Ni)9^as{^^2e,Ni  ®  ^26,^2  )j 

■'®e0[(^2e',V,  ®  *^26,^2  )'^ai(^2e,V,  ®  ^26.^2)] 

“[(‘^2e,W,  ®^2e,N2)9aa{^2e,N,  ®  *^26,^2 )][(^2cV,  ® ‘^2i,V2 )'^aj(^2e,V|  ®  ^2e,V2 )] 

■■*00[(‘5'2i,W,  ®  ^26,^2  )^ra(^2e,V,  ®  ^26,^2)] 

-[fe  ®<^2U)9aa{C2e,N,  ®  •S2e,V2  )][(C2tv.  ®  *S2'.U  ®  ^2..^,)] 

■Ree[[C2U  ®C'2;UK(C2e,V,  ®C2e,V2)]"-  •  • 
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[(^2e,A^j  ®  ®^2e,Ar2)][(  ^2e,N^  ®C2eUK(<^2..A,,®C2..;vJ] 

®^2eU)  "^aa  (^2g,A^2 

“[(^2eVi  ^t7^*^2e,Ni  ®C'2e.A:,)][(  •^2e[Ni  ®<^2;U)  ^sa(^2e,Ni  ®<^2e.wJ] 

'^e^^2e,N^  ®‘y2'«U)  ^as{f'2e,N^ 

■“[(^2e,Wi  ®  *^26,7/2)^^  (*^2e,iVi  ®  ^2e,iV2  )][(‘^2e,^i  ®^2e,A^2)  ^sa{^2e,Ni 
■■®^(?[(‘^2e,A^,  ^sa{f'2e,Ni 

“[(^2e,Wi  ®‘S'2'eUK«(*y2..;..®C'2..;.J][(  *^2e,A^i  ®C^U)^sa(  ^2e,Ny  ®  ^2e,A^2  )] 
‘^0^^2e,N^  ®  ^2e,N^'^ ss{f'2e,N^  ®  ^22,7/2)] 

+[(*^26^1  ®  ^2e,Ni  )‘^aa{f^2e,N,  ®  ^2e,Wj  )][(  ^2e,JVi  ®C2U)  ^ss{f'2e,Ni 
‘^^^[(*^2e,7/j  ®‘y2;U)  ^aa{^^2e,Ni  ®^2e,N,)J 
■^[(‘^2e,7^j  ®S^U)Za  {^2e,Ni  ^2eVi  ^2e,Ni  ®  ^2e,N2  )] 

®s;U)  '^as{f'2e,Ni  ®<^2e.N,)J 

+[(‘^2e,;Vi  ®  *^2e,A^2  )^£ia{f^2e,Ni  ®  ^2e,N2  )][(^2e,;\^i  ®q.UK(  ^2e,7/i 

®QeUK(C2e.2..®C'2..A.,)]" 

‘^[(*^2e,JVi  ®  ‘^2e,Af2  )^a(^2e,W|  ®  ^26,^2  )][(  ^2e,Ni  ®<^2;U)  ®<^2..A^2)] 

®QeUK  (^2e,A^i  ®^2.,iV2)]'' 

*^[(*^26,7^1  ®C^U)Zs[S2e,N,  ® ‘S2e,;.2)][(‘^2tA^,  ® ‘y2;V2K(  ^2e,A7i  ®  ^2e,A^2  )] 
*^^^[(*^2i,7/i  ®‘S'2;U)  ^aa{f'2e,N^ 

®  ^2e,N2  y^as{^2e,N^  ®  *^26,7^2  )][(‘^2e,A^i  ®  ^2e,N2  )^aa(  ^2e,7Vi  ®  ^26,7/2)] 
®  ^2e,N2  )‘^as{^2e,N,  ®  ^26.^2 )] 

“[(*^2e,Af,  ®  ^2e,N2  (*^2e,W,  ®  *^26,^2  )][(‘^2eV|  ®  *^26^2 )  ®<^2e.;^2)] 

■^6>6>[(‘^2e.iVi  ®  ^2e,N2  )^Ja(^2e,7/,  ®  ^2e,N2  )] 

“[(‘^2i,7Vi  ®  ^2e,Af2  )'^as{  *^2e,7V,  ®‘^2e.^2)][(  ^2e,Ni  ®‘y2'eUK«(  ^2e,7^i  ®<^2e.y^2)] 
'^e^^2e,Ni  ®  ^2e,N2  )l^ss  if'2e,Ni  ^^26,^2)]  •  '  • 
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^2e,Ni  ^2e,Ni  ®‘S2’eUK(  ^2e,N^  ®  ^2e,//2)] 

^aa{^^2e,Ni  ®^2e.N^)J 

“[(^2eVi  ^as{f'2e,Ni  ®*2..».)I(C2-,\®-S2'',».)  ^as{^^2e,Ni 

^as{f'2e,N^  ®  ^2e,Ar2)j 

“[(^2e,iVi  ^as{f'2e,Ni  ®  ^2e,N2  )][(^2e,Ari  ®  *^26^2  ^'^as{f'2e,N^  ®  ^2e,Nj  )j 

^sa{^2e,N^ 

“[(^2e,A^i  ®C^U}9asi  ^2e,Ni  ®‘S2.,A,,)][(  ^2e,Ni  ®‘S2;UK(  ^2e,Ny 

®^le,N^  {^2e,N^  ®  ^2e^^  )| 

+[(*^26^  ®-s2;U)  ^as{^2e,Ni  ®  ^2e,N2  )][(‘^2e.Af|  ®  ^2e,N2  )  ^sa{f'2e,Ni 

•^4‘^2;U®‘S2.UK(  ^2e,Ni  ®  ^2e,N2  )] 

+[(*^2e,//i  ®  *^26.^2  )^ai(‘^2e,W,  ®  ^2e,Af2  )][(  *^2e,JVi  ®C2U)  ^sa{f'2e,Ni  ®C2.,A/J] 

■■^5^>[(^2e,W,  ®  ^2e,N2 )  '^as{^^2e,Ni 

+|(‘^2e,Afi  ®  *^2e,//2  y^as{^2e,N^  ®  ^2e,A^2  )][(  ^2e,N^  ®C^U)9^sa{C2e.2^®C,^.f,,)] 

®C’2e.A.2)]" 

+[('^2J,7Vi  ^as(^*^2e,Ni  ®  ^2e,N2  )][(  ^2e,Ni  '^sa{f'2e,Ni  ®<^2.,Ar2)] 

*^<96>[(^2e,A^i  ®C^Uys(C2e,N®C,^,.,)f 

®‘S'2;U)  ^as{f'2e,N^  ®  ^2e,A^2  )][(  ^2e,JVj  ®  ^2e,Ni  )] 

'^aa{f'2e,N^  ®^2..Af:)]^ 

{f'2eMi  ^2e,N^  ®C'2;UK(<^2e.A..®<^2e.i.2)] 

'^G^f'2e,Nx  {f'2e,Ny 

+[(^2e,Afi  ®  *^26^2  )  ^as{f'2e,Ny  ®C2e,N,  )][(^2e,^i  ^ss{^2e,Ni  ®C'2..W2)] 

’^90^^2e,Ny  ®C2U)9‘sa(C2e,N®C2e,N,)J 
+[(^2eVi  ®-s2;U)  ^as{f^2e,Ni  ^2e,Ni  ®  ^2e,N2  )  ^ss{f'2e,Ny 

■■®e^[(^2e,Wi  ®C2U)9^ss(  ^2e,Ni  ®C2e.N,)f--- 
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[(^2;^ 

•^..[(*52;^  ^S^U)^4^2e.N,  ®C'2.,;.J]" 

-[(C2U  ^S2U)M^2,.^  ®‘y2e.;.J][(^2’..  ®^2-^.;v.K(c,,^, 
■Roe[(C2U  ®Su,N,yas[C2.,N, 

-[(^2^.,  ®^2-J...K(*y2...  ®-y2...)][(*s2;'...  ®^2;UK(c2,.. 

■Re0[(S2U  ®C'2;UK(<^2..Ar, 

~[(C2U  ®‘y2..J][(^2.U  ®‘S2t..K(C2,;.. 

•Ro0[(C2U  ®C^U)M^2e,N, 

+[(*^26,^,  ®  ^2e,N^9sa{C2e.Ni  ®  *^2e,yVj  )][(^^2i.Af,  ®  ^2l,N^9‘„^{C2^  ®C’2e.Afj)] 

•^..[(‘y2tA.,  ® ‘y2t;..K(C2.,A,, 

+[(</.,  ®  *S'2;UM.(C2.,^.  ®*52..A/J][(C2-e'A,.  ®  <52;^  K(q,,;V.  ®C'2..A.J] 

•^.^(^2;^.  ® ‘S'2;UK(c2,,^. 

+[(‘S2;'a,,  ® ‘S2;^,)?.„(C2,.^,  ®‘y2a.;.J][(C2;U  ®  S2U)9^as{C2e,N, 

■•*5^[(‘5'2eV,  ®^2j.wJ^ia(C2e,Ar,  ®C'2e,JV2)] 

+[(‘^2e,A^,  ®*y2e,ArJ?M(C2e.//,  ®‘S'2^  ;^,J|(C2jj^^  ®  *^2 ®<^2e.W2)] 
®C2U)9‘ss(C2e.N,  ®C2e,A.J]" 

-[(<^2.V.  ®C2e,N,)9a[S2e,N,  ®  C'2.,JV,  )J[(‘S'2t^,  ®  ^2^^  (^2,,^,  ®  <^2..yvJ] 

'^0^^2l,N^  ®‘^2j.Afj)^oa(C2e_Ar,  <8) 

-[(^2^,;/,  ®<^2'.U)?.«(*y2e.^,  ®^2..AfJ][(‘S'2;',w,  ®  <^2'.^  K«(<^2e,VV,  ®  C2,_;^Jj 

®  •S2;UK(C'2e.;V.  ®C2.,N,)f 

~^2e,N^  ®^2e,Wj)?ra(‘S'2e,Af,  ®C'2g  ;vJ|(52j_^^  ®^2l,N^'^sa{C2e,N^  ®C2c,yVj)] 

■-*^^[(*^2^.^?,  ®^2l,Ny^sa{C2e.N,  ®C2e,N^ 

~]^2e,Ny  ®^2e,N^)9sa{S2e,N^  ®  ^2e,AA2  )][(‘S'2  J,Af,  ®  ^2l,Ni)'^sa{f2e,N,  ®^2e.A^2)] 
'^e^2l,N^  ®  ^2eU2  )‘^!:s{f2e,N^  ®  ^2e,N^  )]  ’  '  ' 
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{^2e,Ni  ^2e,Ni  ^2e,Ni 

'^Oe^{^2e,Ni  '^aa{f'2e,Ni  ®  ^2e,N2  )] 

®  ^2e,N2  )^sa{^2e,Ni  ®  ^2e,N2  )][(  ^2e,Ni  ^ss{^^2e,Ni  ®  ^26,^2)] 

'-®06>[(^2e,A^i  ® ‘^2eV2)^«*y(^2e,A^i  ®^2e,//2)] 

■^[(‘^2i,iVi  ®C^U)9sa(  ^2e,Ni  ^2e,Ni 

'^a0^{^2e,Ni  ^saif'le^Ni  ®C2e,N,)J 

"•■|(*^2e,iVi  ®C^U)9sa  {f'2e,Ni  ®  ^2e,N2  )][(^2e,i\^i  ®C^U)9^ss(  ^2e,Ni  ®  ^2e,N2  )] 
®  ^2e,N2  )^ii(^2e,^/i  ®  ^2e,H2  )j 

+[(‘^2i,//,  ®  ^2e,N^  )5ii(‘^2e,A'|  ® ‘^2e,wJ][(‘^2i.W,  ®  ^2l,N^9faa{f'2e,N^  ®  ^26,^2)] 

■■®6'6'[('^2e,W,  ®  ^le,N2  )^aa(^2e,//i  ®  ^2e,/^2  )] 

■^[(‘^2e,//i  ®*S2tiV2)  ®  *^26,^2  )][( 

'^e^^2e,N^  ®‘S'2;V2)  ^as{^2e,N^  ®<^2e,Af2)]'^ 

+[(*^2e,W,  ®  ^2e,N2  (*^2e,A^i  ® ‘^2e,//2)][(*^2e,A^i  ®  ‘^2e,/^2  )  "^^70  (^2e,A^i  ®f^2..A.2)] 

®c2;UK  (^2e,A^i 

+[(‘^2e,Afi  ^ssi^2e,Nx  ®‘^2.,;v2)][(  ^2e,N^  ®  ‘^2cV2  )  ^aa{^^2e,Ni  ®  ^2e,N2  )] 

'^e^^2e,Ni  ®  ^2e,N2  ^'^ssif'le.N^  ®  ^2e,N2  )j 

■*"[(^26,^^!  ®'S2''.».)  ®  ^2e,N2  )][(^2e,A^i  ®  ^2e,N2  )'^as{^2e,Ni  ®  ^2e,N2 )] 

®  ^2e,N2  )^aa  (^2e,Ni  ®  ^2e,N2  )] 

®<7.2K(<^2.,;.,®'52..2V2)][(  ^2e,Ni  ®^2-eU)  ^as(^2e,Ni  ®<^2..AA2)] 
®S2lN^)  ^as{^2e^Ni 

+[(^2e,Afi  ®S2U)  ^ss{f'2e,N^  ®  ^2e,N^  )][(^2e,W|  ®  *^26,^2  ) 

■^^^[(*^2e,A^i  ®  ^2e,A^2  y^sa  {f'le,^^  ®  ^2c,7\^2  )] 

■‘"[(^2e,//i  ®‘S2;U)  ^ssiple^N^  ®  ^2e,N^  )][(^2i,W,  ®  ‘^28.^2  )^oi  (  ^2e,N^  ®C'2.,y.2)] 
■■*S^|(^2e,W|  ®C2U)  ^ss[^2e,Ni  ®C2e,N,)f--- 
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^2e,Ni  '^sa{f'2e,N^  ®<^2e,A^,)] 

■•*es[(‘^2e,W,  ^aa{f'2e,N^ 

+[(‘^2e,Ar,  ®C^U)9ss(  ^2e,Ni  ^2e,Ni 

'^e^^2e,Ny  ^2e,Ni 

+[(*^2e,Af,  ®<^2’eVj  ^ss{^2e,Ni  '^sa{f'2e,Ni 

*^^^[^(*^2e,//i  ®  ^2e,N2^‘^sa{^2e,Ni  ®^2e,;'/2)j 
+[(‘^2e,//i  ^^2e,N2)^ss{^2e,Ni  ®^2e,i^2)][(  ^2e,Ni  ^sa{f'2e,Ni  ®<^2e.W2)] 

■■^^^[(^2c,W|  ®<^2;U)  '^ss{f'2e2N^ 

+[(^2e,//i  ^ss{f'2e,N^  ^2e,iVi  ®  ^2e2N^'^ ss{f'2e2N^  ®  ^2e,A^2  )J 

®  ^2e,N^  y^aa  {f'2e,Ny  ®  ^2e,Ar2  )  j 

+[(^2e,//i  ^ss{f'2e,Ni  ®<^2..W2)][(  ^2e,Ari  ®C2;U)  '^ss{f^2e,Ni  ®  ^2e.N2  )] 

■•®w[(^2e,JV,  ®  ^2l,N2)'^as{f-2e,Ni  ®^2e,Af2)] 

■•■[(^26,^1  ®^2l.N2)9ss{^2e,Ni  ®  ^26,^2  )][(^2i,JV,  ®  ^2e,Af2  )^ii(^2e.;\/,  ®  ^26.^2  )] 
■•®^^[(‘^2e.Af,  ®  ^2e,N2  )  '^sa{^^2e,Ni  ®^2e,N2)f 
+[(^2e,Af,  ®<^2;U)  *^ss{^2e2N^  ®  ^2e,A^2  )][(  ^2e,iVi  ®C^lN2)‘^ss[C2e,2^®C^e,N2)] 

®  ^2e,/^2  )'^ss{f'2e,Ni  ®  ^2e,N2  )] 

+[(^2eVi  ®C2-.V.)?*(  ‘^2e,A^i  ®  *^2e,A^2)]  ®  ^2e,N2  )  (‘^2e,Af,  ®  *^26,^2  ) 

“[(^2e,Af,  ®C^U)9as  (^2e,A^i  ®  ^2e,N2  )j^w  {f'2e,Ni  ®  *^2e,A^2  )  {f'2e,N^  ®  *^2e,//2  ) 
’“[(^2e,^i  ®  ^2e,;\^2  ®  ^2e,A^2  )]-^w(*^2e,A^i  ®  ^2e,A^2  )  (*^2e,;\^,  ®  ^2e,A^2  ) 

+[(^2e,A^j  ?w(^2e,//i  w  (^2e,^j  ®  ^26,^2)  (^2e,Ar,  ®  ^26,^2  ) 

“  if'2e,Ni  ®  ^2e,A?2 )  {f'2e,Ni  ®  ^2e,N2  )-®^^[(  ^2e,Ni  ®  ^2e,N2^‘^aa{^2e,Ni  ®  ^26,^2)] 
+(^2e,W,  ®  ^2e,N2  )  {^2e,N,  ®  ^2e,N2  )-®«2[(  ^2e,Ni  ®  *^2], ^2  y^asi^  ^2e,Ni  ®  ^2e,N2  )] 
+(^2e,W,  ®  ^2e^2  )  (^2e.W|  ®  ^2e,N2  )'®^«'[(  ‘^2e,iV,  ®C2;U)  ^sa{p2e,N^  ®^2e.N2)( 
'^{f'2e,Ni  ®  ^2e,N2 )  (^2e.Af,  ®  ^2e,Af2  )-®6’^[(  ^2l,N^  ®C^In2)  ®<^2..Af2)]'"- 
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Equation  (156)  can  be  simplified  by  making  the  substitutions 


®  ^2e,N2  ®  ^2e,A^2 )  ’  (1  ^^7) 

^{^2e,Ni  ®  ^2e,N2)^ww{^2e,N^  ®  *^2e,iV2  )  ’  (1^^) 

^'^as^as  ~  {^2e,Ni  ®  ^2e,N2  ^^ww{f'2e,N^  ®  *^2e,A^2  )  ’  ( ^  ^9) 

^'i(^sa'^sa~{^^'2^^^^i®^2,e,N^^ww{^^2e,N^®^2e,N^  ?  (160) 

”  {f'2e,N^  ®  ^2e,N2  ")^ww{f'2e,N^  ®  ^2e,A^2  )  *  (161) 

Making  the  substitutions  from  Eqs.  (157)  -  (161)  into  Eq.  (156)  produces 


{^2e,Ni  ®  ^2e,N2  )  ^aa'^aa^0^@^,^aa  {^2e,Ni  ®  ^2e,N2  ) 

'^{f'2e,Ni  ®  ^2e,N2  )  ^aa^aa^0,,e,,'^as{f'2e,N^  ®  ^2e,N2  ) 
'^{f'2e,N^  ®C2e,N,y9aa‘^aaR.  9,fi,,'^sa{^2e,N, 

®  ^2e,Ni )  ^aa^aa^e^0^^ss{^^2e,N,  ®  ^2e,N^ ) 

“(‘^2e.A',  ®  ^2e,N2  )  '^aa^as^0,,0,,'^aa  ('^2e,A^i  ®  *^2e,A^2  ) 

“(*^2e,JVi  ®  ^2e,N^  )  5^aa^oi'^0„©„^oi(^2e,W,  ®  *^26,^2 ) 

“(*^26,^1  ®^2e,A^2) 

“(*^2e,A^i  ®  ^2e,A^2  )  as^0,(0,^ !is{^2e,}^^  ®  ^2e./^2  ) 

~(f^2e,Ni  ®  ^2e,Ni )  ‘^aa'^sa^@^@,,'^aa{^2e,N^  ®  ^2e,N2  ) 

~{^2e.Ni  ®  ^2e,N^  )  ‘Paa^sa^0,,@,,‘^as{f'2e,Ni  ®  ^2e,Ni ) 

~{f'2e,N^  ®^2e,N2)  "^aa^ sa^0J9,^sa\^2eMi  ®^2e,N^ 

~{p2e,Ni  ®  ^2e,Ni )  %a^sa^0^0j^ss{f'2e,Ni  ®  ^2e,N2 ) 
+(*^26, Af,  ®  ^2e,Ni  )  ?oa^ji'®0„0„^aa(‘^2e,W,  ®  ^2e.N2  ) 

+(*^26, Wi  ®‘^2e,A'2)  7aa^ss^0^0,,^as{^2e,Ni  ®S2e,N2y'^ 

'^{^2e,Ni  ®  ^2e,N2 )  ?aa^Ji'®0„0„'^ia(‘^2e.//|  ®  ^2e,N2  ) 

+(‘^2e,2/i  ®  *^26,^2)  %a^ss^0,,0J^ss{f'2e,N^®^2e,N^  • 


no 


®  ^2e,W2  ^ 

)  %s'^aa^0^,.e,^^aa  (*^26, AT,  ®  ^2e,N^  ) 

/ 

i-1 

/  \-r 

~\^2e,Ni 

/  'Pas^aa^0^e^^'^as\f'le2N^  ®  *^26,^2  j 

/ 

i-1 

/  \-r 

-(^2e,N, 

/  'Pas'^aa^0,,0,,'^ sa\^2e,N^  ®  ^2e,//2  j 

/ 

1-1 

/  x-r 

®<^2..Ar2’ 

/  ss\f'2e,N]^  ®  ^2€,JV2  j 

v-1 

/  \-r 

~{^2e,Ni 

®C’2e,W2' 

)  ^as^as^0^J0^^aa  (‘^2e,JVi  ®  *^2e,JV2  j 

( 

v-1 

/  x-r 

~\^2e,N, 

®t’2.,2/2' 

)  ^as'^as^0^J0^^as\p2e2Ny^  ®  *^2e,A^2  j 

/ 

v-1 

/  x-r 

\^2e,Ni 

)  ^a^as^0J0^^sa\^2e,Ny^  ®  ^2e,A^2  j 

/ 

v-1 

/  x-r 

\^2e,Ni 

®C2e,W2' 

)  as^0J0^^ ss\^2e;M^  ®  ^2e,Ar2  j 

/ 

.-1 

/  x-r 

+(*^26,//, 

®‘^2.,W2; 

1  %s^sa^0,J0,^aa  (*^2^,//!  ®  ^2e,N^  ) 

®‘52e,iV2; 

1  ?fli^sa'*0„0„^a4^2e.W,  ®*^2e,Wj) 

/ 

.-1 

/  x-r 

■'■{‘^28, W, 

1  9as^sa^0^^0,,^sa  (*^2e,A^i  ®  Qe,A^2  j 

/ 

.-1 

/  x-r 

+(‘^2e,W, 

®‘y2.,;.2; 

I  ®  ^2£;,//2  j 

■•■(^26,2^1 

®‘^2..2V2' 

)  ^as*^ss^0^,0,,'^aa  {^2e,Ni  ®  ^2e,N^  ) 

/ 

> 

,_1 

/  x~r 

+(^2e,W, 

®‘^2e.A22, 

)  ‘?as'^ss^0,,0^,'^as\f'2e,Ni  ®  j 

/ 

> 

i-1 

/  x-r 

)  ^as^ss^0^^0,,^sa  (*^2e,A^i  ®  ^2e,N2  ) 

/ 

> 

.-1 

/  x-r 

)  \^2e,iVi  ®  ^2e,A^2  j 

> 

,„1 

/  x-r 

“(^2e,A2, 

®*52..W2, 

)  7sa'^aa^0,,0,^^aa\^2e,N^  ®  '^2e,A^2  j 

/ 

> 

.-1 

/  x-r 

“(^2<;,A2, 

)  ^sa^aa^0J0,^as\^2e,N^  ®  *^2e,A^2  j 

“(^2e,A2, 

®^2e,N,] 

)  ?ra^aa-®0^,0„^ia(*^2e,A/,  ®^2e,N^ 

\ 

/  x-r 

®^2e,N,, 

)  '?sa'^aa^0^,0,^^ss\f'2e2N^  ®  ^2e,//2  j 

+(‘^2e,W, 

®^2e,N,] 

1  ysa‘^as^0^0,,‘^aa{^2e,Ni  ®  ^2e,N^ 

+(‘^2e.Wl 

*^^2e,N^] 

1  7sa‘^as^0^0j^as{f'2e,N2  ®  ^2e,N2  ) 

/ 

> 

.-1 

/  x-r 

®^2e,N,^ 

I  5^sa^a^^0„0„^^a(‘^2e,A^i  ®  ^2e,N2  ) 

/ 

> 

.-1 

/  x-r 

■'■(‘^26, W, 

®'^2..W2, 

I  ^sa'^as^0,,0,,^ss\^2e,N^  ®  ^2e.7^2  j 

/ 

v-1 

/  x-r 

“(^2e,W, 

®  ^2e,N2 

)  9sa‘^sa^0^0,,  ^aa  (*^2e,A^i  ®  ‘^2e,A^2  / 

/ 

v-1 

/  x-r 

“"(^2e,JVi 

)  ®  ^2e,JV2  j 

Ill 


~{f^2e,Ni  ®  ^2e,Arj  )  ^sa'^sa^0^0^‘^ssif'2e,N^  ®  ^2e,N2  ) 
+(‘^2e,Af,  ^^26,^2)  ^sa‘^ss^0^0^^aa{^2e,Nt  ^•^26,7/2) 

+('^2e,Af,  ®  ^2e,Ni )  ^sa^ss^0„0j^as{f'2e,Ni  ®  *^2e.JV2 ) 
■'■(‘^2e,Afi  ® ^26,^2 )  ^sa^ss^0^0j^sa{p2e,Ni  ®  ^2e,A72 ) 
+(‘^2e.W|  ®  ^2e,N^  )  ‘Psa‘^ss^0„0,,  '^ss(^2e,N,  ®  ^2e,N2 ) 

+(*^2e,Af,  ®  *^2e.Af2 )  5^ii^aa'®6)„0„'^aa(*^2e,//,  ®  ^2e,N2  ) 

■'■(‘^2e,Afi  ®  *^26.7^2  )  ^ss^aa  ^0^0,,‘^as(^2e,N,  ®^2e,2^2)'^ 
+(*^26,7/,  ®  '^2e,N2 )  ^ss'^aa^0„0,,^sa{^2e,Ni  ®  ^2e,N2 ) 

■'■(‘^26,7/,  ®  ‘^2«.A'2  )  ^ss'^aa^t  ^^0j^ss{f^2e,Ny  ®^2e,N^ 
"^(^2e,JVi  ®  ^2e,N2 )  '^ss^as^0,,0^,^aa{^2e,Ni  ®  *^2e,A^2  ) 

~^{f'2e,Ni  ®  ^2e,N2  )  '^ss^as^0,,0,,^as{f'2e,Ni  ®  *^2e,^2  ) 

'^{f'2e,N^  ®  *^2e,//2  )  ^ss^as^0^0^^'^sa  {^2e,Ni  ®  ^2e,N2  ) 

~^{^2e,Ni  ®  ^2e,N2  )  ^ss^as^i  0„0„^«(<^2e.7V, 

~^{^2e,Ni  ®  ^2e,N2  )  ^ss^sa^0^0^, ^aa (*^2e,//i  ®  ^2e,N2  ) 

‘^(*^2e,A^i  ®  ^2e,N2  )  ^ss^sa  ^@^@,,‘^as{f2e,N,  ®S^e,N2Y 
■’■(‘^26,7/,  ®  ^2e,N2 )  “^ss^sa^i  9^0j^sa[S2e,N^  ®C2e,N2Y 

“*'(‘^2e,JVi  ®C’2..7V2r'5'A^^  i^0,,‘^ss^2e,N^  ®^2e,N2Y 

+(^2e,7V',  ®  ^2e,N2  )  ^ss^ss^i  9^0,,^aa{^2e,Ni  ®S2e,N2Y 
+(^2e.W,  ®  ^2e,N2 )  ^ss^ss^i  9^,0,,!^as{f'2e,N^  ®^2e,N^ 
■*■(^26,7^,  ®  ^2e,N2  )  ^ss^ss^i  9^^0^,^sa{^2e,Ni  ®^2e.N2Y 
+(^2e,W,  ®  ^2e,N2 )  ^ss^ss  ^0^0,,‘^ssY'ie,N2  ®^2c,nY'^ 
■*■(^26,77,  ®  ^2e,N2  )  ®  *^26.7^72) 

~(^2e.W,  ®  ^2e,A^2  )  ^  *^26,7^2  ) 

“(^2e,^i  ®  ^2e,^2  )  ^sa^TOsa'^sa  (*^2e,A^i  ®  ^2e,N2 ) 

■'■(^26.77,  ®  ^26,7^2)  ^s.s'®»t,»tj(^2e,77,  ®  ^6,772)  ■•• 
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(162) 


®  ^2e,N-i )  ^0,,@J^aa  (‘^2e.W,  ®  ) 

+(^2c.W,  ®^2e,Nj  ^&^@J^as{f2e,N^ 

'^{f'2e,N^  ®^2e,N^  •®(9j,(9„^ra(‘^2e,W,  ®  ^26,2^2) 

+(^2£!,Afi  ®^2e,N^  ^@,,0j^ss{f'2e,Ni®^2e,N^  ■ 


The  simplification  of  Eq.  (156)  to  produce  Eq.  (162)  makes  use  of  the  facts  that  =  r?(^. 


‘^as  =  ‘^as^  '^sa=‘^sa^  and  1^ss  =  ‘^ss  because  the  matriccs  and  are  all  di¬ 

agonal.  Combining  terms  in  Eq.  (162)  results  in  the  expression 

{f'2e,N^  ®  ^2e,N^ )  ^aa^aa^B^eJ^aa  ~  5oi^oj-®0„0„^oa  ■•■  ^aa^Ki„^7ci^„ 

“5^ia^ra'®0„0,,^aa  +  ?M'^ii'®0„0„'^aa](‘^2e.W,  ® ‘^2e,2V2) 

+(^2e,W|  ®  ^2e,N2  )  [?aa^aa-*0„0„^aj  “  ^as^as^B^eJ^as  “ 

~?ia'^M-®0„0„^ai  +?M^M'®0„0„^ai](^2e.W,  ®  *^26,^2) 

+(^2<;,2V,  ®  ^2e,N2  )  [paa'^aa^B^B^^sa  ~  Pas'^as^B^B.^^sa  ~  Psa^W^TV^ 

~Psa‘^sa^B^B^^sa'^PssP^ss^B^B„P‘sa^{^2e,Ni  ®  ^28,^2) 

+(^2e.Afi  ®  ^2e,N2  )  [5^aa^aa'®0^0„^ii  “  PosP^qs^B^bJ^ss  +  Pss^K/JO^ 

~PsaP^sa^B^JP^ss'^PssP^ss^B^B,P^s^^2e,N^  ®^2e,N^ 

+(*^26,^,  ®  ^2e,W2  )  [“5^aa^ai-®0„0„^aa  ~  PasP^aa^B^B,P^aa 

'^Psa^ss^BjBj^aa'^Pss'^sa^B,,B,P^aa^^2e,Ni  ®  ^2e,N^ 

+(‘^2e,W,  ®  ^2e,N2 )  [“5^aa^aj-^0„0„^ai  “  PasP^aa^B^B,P^as 

'^Psa^ss^B^Bj^as  +  ?si^ia'*0„0„'^a^](^2e,W,  ®  ^2e,N2 ) 

+(*^2e,#,  ®  ^2e,N2  )  [~5^aa^ai-®0„0„^ia  “  Pas'^aa^B^B.P^sa 

+?ra^ii'®0^0„^ia  +5ii^sa'®0„0„^ia](‘^2e,2V,  ®  ^26,^2) 

+(52,,;,,  ®C^e,N)\9aaP^asRB^Bj»ss-PasP^aaRB^B^P^^^ 

-^PsaP^ss^B^Bj^ss'^PssP^sa^B^B,P^s^2e,Ny  ®^2e,N^ 

+(^2e,W,  ®  ^2e,N2  )  \~PaaP^sa^B^Bj^aa  +  PasP^ ss^B^bJ^oo 

"PsaP^aa^B^B.P^aa  +  ?ii^as^0„0„‘^aa](*^2e,#i  ®  ^2e,N^  •  '  ' 


\-T 


\-T 


\-T 


\-T 


\-T 


\-T 


\-T 
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+(^2e,JV,  ® 


+1 


’26,^2)  [  ^aa'^sa^e^ej^as'^^as^ss^e^ej^as 

-^sa‘^aa^0^0,,^as  +  ^ss'^as  •®0„0„^as](C2e,W,  ®‘y2e,Af2)"'' 
\^0^0j^sa  +  ^as^ss^0„0,,^i 


^  sa  -  aa  -  as  ''ss^as  -  -  as ze,j\^ 

'{ple,Ny  ®  ^le,N^ )  \^^aa^sa^0^0j^sa  +  ^as^ss^0„0,,^sa 

~^sa‘^aa^0^0,,^sa  +  ^ss'^as^0,,0,,^sa^^2e,Ni 

/  \-l  r 


-‘Psa‘^aa^0^0,,‘^sa'^9ss'^as  ^@,s0ss^sa\{^2e,Ni  ®C2e.MS 
'^{f'2e,N^  ®^2e,N^  [  'Paa^sa^0^0^ss'‘^'?aP^ss^0jaJ^ss 

-‘Psa'^aa^0^0j^ss  +  ^ss'^as  ■*0„0,,^5i](C2e,W,  ®C2e,Ny 

+('^2e,A^i  ®  ‘^2e,Ar2  )  ^aa‘^ss^0^J^aa  +  ‘Pas‘^sa^0^0^‘^aa 

+?ra^oi-^0„0j,^aa  +?«^ao'^0^®,,^aa](‘^2e,W,  ®  *^26,^2) 
+(*^26,^,  ®  ^2e,Ni  )  [9aa'^ss^0^0„‘^as  +  ^as‘^sa^0^0^^as 


+7sa'i^asJ^0^0j^aa+9ss’^aaJi0^0,:’*aal^2e,N,  ^2e,N,) 

^(*^26,^1  ®  *^26.^2  )  \_‘^aa'^ss^0^0„^as  +  ^as^ sa^0^0j^as 

'^'?sa^as^0J0,^as  ^ss^aa^0^0^^asJ(^2e,N,  ^  ^2e,N2^ 

^{^2e,Ni  ®  ^2e,N2 )  ^aa^ss^0^0j^sa  +  ^as^sa^0^0j^sa 

+5^sfl^ai-®0j,0j,^sa  +?ii^aa^0„0„^io](‘^2e,W,  ®  ^26,^2) 
+(‘^2e.Wi  ®  *^22.^2  )  +  9as'^sa^0^J^ss 

r'^0„0„.^ii  +  ?M^ia-^0„0,„‘^M](^2e,Af|  ®^2e,A^2) 

®-S2e.;v2)'" 


-T 


-T 


-T 


'{f'2e,N^®^2e,N^ 

(^2e.W,  ®  ^2e,N^ ) 


j,N^)  ^0^0^‘^aa{^2e,N^ 

+(^2e,W,  ®  Qe.W2  )  ■^0„0„^ot(^- 
+(^2e,W,  ®  ^26,^2 )  ^0^0„^sa{^2e.Ni  ®  ^2e,Ar2 ) 
+(^2e,Wi  ®C’2e.iV2r^,  0„0„^«(^2e,Afi  ®  Qe,A^2  ) 


-T 

T 

-T 

-T 


Equality  in  Eq.  (163)  will  hold  only  if 


^aa{^aa^e,,0s,'^aa'^  ^as^as^0^^©^^aa  ^sa^sa^0^,0,,'^aa 

'^^ss*^ss^0^^0ss^aa  "*  ^6>j,0j_y^aa  ’ 

^aa^aa^0^jB^^ as  ~  as^0 ,J0  as  )  “  ^sa^sa^0^^0^^'^as 

'^7ss*^ss^0,,0j^as  ~  ^0,,&^^asy 

'?aa^aa^0,<0,!^sa  ~  9as^ as^0^J0,^ sa  ~  "9 sa{^ sa^0 ,Jd sa  +  ) 

~^"Pss*^ ss^0^,0,!^ sa  -  ^©..©..'^sa^ 


(163) 


(164) 


(165) 


(166) 
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(167) 


‘9aa‘^aa^0,,@J^ ss  %s‘^as^@,,0j^ss  7sa‘^sa^0,,0j^ss 

'^^ss{^ss^0„0,,‘^ss  ■*■  ^7C^7li,)~  ^0^0,,^ss’ 

~^aa^as^0^0j^aa  ~%s^aa^0^0j^aa'^^sa^ss^0^0j^ac 
^^ss'^sa^0^0j^aa  ~  0> 

-‘Paa‘^as^0,,0j^as  “  ^cJ^aa^0^0„^as  +  9sa‘^ss^0,.0^^as 
'^‘Pss^sa^0,,0j^as  ~  ®> 

~%a^as^0,,0j^sa  ~^as^aa^0^0^^sa  ^sa‘^ss^0^0„^sa 

+^ss'^sa^0^0^^sa  -  0’ 

~‘^aa^as^0^0j^ss  ~  ^as'^aa^0^0„^ss  +  %a^ss^0^0j^ss 
+^ss^sa^0^0j^ss  =  0’ 

-9aa‘^sa^0,^^‘^aa+9as^ss^0^0,,‘^aa  “  ?sa^<7a-*0„0^^a. 
+9ss‘^as^0,,0j^aa  =  0’ 

~9aa‘^sa^0,,0j^as  +  9as‘^ss^0„0j^as  ~  909^00^0^0  9^ as 

^9ss9fas^0^0j9^as  =  0, 

~9aa9^sa^@,,0j9fsa-^9as9^ss^0^0  9fsa  ~9sa9‘aa^0^0„9^sa 
+9s9‘osR0,,0„9^sa  =  ^, 

~9aa9^sa^0^0j9^ss  9as9^ss^0„0„9^ss  ~9sa9^aa^0^0j9^ss 
+9s9‘osJi0A^ss=^^ 

9aa9^ss^0 ^0^9^00  ■*"  9as9^sa^0,^0j9^aa  "*"  9sa9^as^0„0„9^aa 
+‘Pss9‘aa^0^0^9‘gg  - 

9aa9^ss^0,,0,9^as  9as9^sa^0,^0j9^as  9 sc^ as^0 ^0 ^ as 

'^9ss9^aa^0,,0j9fas  =  0’ 

9oa9fssR0^0j^so  +  9os9‘so^0A9'sa  +  9so9^osR0^0j^ so 

'^9ss9‘aa^@a9‘sa  =  0’ 


(168) 

(169) 

(170) 

(171) 

(172) 

(173) 

(174) 

(175) 

(176) 

(177) 

(178) 
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and 


(179) 


9aa‘^ss^0^0^^ss  +  ‘Pas‘^sa^0,,0j^ss  +  ^sa‘^as^0^0j^si 

'^^ss^aa^0„0,,^ss  ”  ®- 


Just  as  in  the  one-dimensional  derivation  of  the  scalar  Wiener  filter,  every  matrix  in 
Eqs.  (164)  -  (179)  is  either  diagonal  or  approximately  diagonal.  The  matrices  “Psa^ 

9ss’  “^aa’  ^as’  ^sa’  and will  always  be  diagonal.  The  transform  domain  noise  correlation 
matrices,  will  be  diagonal  because  the  noise  samples  are 

uncorrelated  and  have  uniform  variance.  The  matrices  ^  ^  ^  ,  and  R^  ^  are 

*^aa*^aa  *^as*^as  *^sa*^5a  *^ss*^ss 


thus  exactly  represented  by  their  diagonal  elements  70^^{k^,k2), 

The  transform-domain  object  correlation  matrix,  Re^e^,  will  be  well-approximated 

by  its  diagonal  elements  based  on  the  same  assumption  of  a  highly-correlated  object  used  for  the 
one-dimensional  case.  The  matrix  Re^e^  is  therefore  well-approximated  by  its  diagonal  ele¬ 


ments  ©]Xk^,k2). 

The  above  approximations  allow  for  a  scalar  solution  based  solely  on  the  diagonal  ele¬ 
ments  of  the  matrices  in  Eqs.  (164)  -  (179).  Solving  for  the  diagonal  elements  of  Eqs.  (164)  - 
(179)  produces  the  16x4  overdetermined  matrix  equation  which  appears  as  Eq.  (180)  on  the 
following  page.  Each  entry  in  Eq.  (180)  is  dependent  on  the  transform  domain  index  variables 
Atj  and  k^,  but  this  dependence  is  suppressed  in  the  equation  to  conserve  space.  The  16x4  ma¬ 
trix  equation  in  Eq.  (180)  reduces  to  the  7x4  matrix  equation  of  Eq.  (132)  from  Section  4.3.2 
which  is  repeated  here  following  Eq.  (180). 
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increased  savings  in  computational  complexity  for  symmetric  point  spread  functions.  The  fact  that  the  discrete  Fourier 
transform  diagonalizes  a  circulant  matrix  provides  an  alternate  way  to  derive  the  symmetric  convolution-multiplication 
property  for  discrete  trigonometric  transforms.  Derived  in  this  manner,  the  sym-metric  convolution-multiplication  property 
extends  easily  to  multiple  dimensions  and  generalizes  to  multidimensional  asymmetric  sequences.  The  symmetric 
convolution-multiplication  property  allows  for  linear  filtering  of  degraded  images  via  point-by-point  multiplication  in  the 
transform  domain  of  trigonometric  transforms.  Specifically  in  the  transform  domain  of  a  type-II  discrete  cosine  transform, 
there  is  an  asymptotically  optimum  energy  compaction  about  the  low-frequency  indices  of  highly  correlated  images  which  has 
advantages  in  reconstructing  images  with  high-frequency  noise.  The  symmetric  convolution-multiplication  property  allows 
for  well-approximated  scalar  representations  in  the  trigonometric  transform  domain  for  linear  reconstruc-tion  filters  such  as 
the  Wiener  filter.  An  analysis  of  the  scalar  Wiener  filter's  improved  mean-squared  error  performance  in  the  trigonometric 
transform  domain  is  given. 
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